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1  Introduction 


Chemical  lasers  represent  an  important  component  in  many  advanced  weapons  systems  of  interest  to 
the  Air  Force,  e.g.,  the  airborne  laser  program  (ABL).  The  chemical  oxygen-iodine  laser,  or  COIL, 
is  a  short  wavelength,  high-power  chemical  laser  that  operates  on  an  atomic  iodine  laser  transition. 
COIL  was  invented  in  1977  at  the  Air  Force  Weapons  Laboratory,  now  a  part  of  the  Air  Force 
Research  Laboratory  at  Kirtland  Air  Force  Base,  New  Mexico.  AFRL/DE  continues  to  develop 
COIL  technology  using  as  a  primary  test  apparatus  the  10  kilowatt-class  device  called  RADICL 
(Research  and  Development  Iodine  Chemical  Laser) 

Research  and  development  programs  for  COIL  typically  employ  numerical  analysis  techniques 
such  as  computational  fluid  dynamics  (CFD)  to  provide  insight  into  laser  performance  and  design. 
A  common  goal  shared  by  these  programs  is  the  need  to  more  thoroughly  understand  the  underlying 
physical  processes  of  chemical  lasers  so  that  performance  may  be  optimized.  Enhanced  CFD  analysis 
tools  can  play  an  essential  role  in  this  regard.  High-resolution  CFD  simulations  can  aid  in  interpreting 
and  explaining  measured  data,  make  predictions  for  operating  conditions  beyond  the  test  database, 
and  provide  accurate  estimates  of  laser  power  output.  Beyond  predicting  the  performance  of  a  given 
laser  configuration,  new  design  tools  are  required  to  aid  scientists  in  tuning,  scaling  and  shaping  a 
practical  system.  AeroSoft’s  research  and  development  program  has  focused  on  the  following  three 
goals: 

1.  To  develop  accurate  and  efficient  CFD  analysis  tools  for  modeling  COIL  lasers. 

State-of-the-art  CFD  methods  have  been  extended  to  include  new  capabilities  and  physical 
models  important  in  COIL-type  laser  systems  such  as  the  RADICL. 

2.  To  develop  high-level  software  tools  tailored  for  the  design  and  optimization  of 
COIL  lasers. 

The  continuous  sensitivity  equation  method  (CSEM)  for  CFD-based  design  has  been  extended 
to  include  new  capabilities  and  physical  models  important  in  COIL-type  laser  systems. 

3.  To  develop  a  framework  for  coupling  discipline-specific  sensitivity  analyses. 

Work  specific  to  this  contract  has  led  to  new  sensitivity-analysis  tools  for  the  coupled  fluid- 
dynamic/power-extraction  systems  in  COIL-laser  designs. 

The  Phase-II  COIL  effort  has  focused  on  extending  AeroSoft’s  GASP  and  SENSE  software  for 
COIL  lasers.  GASP  is  a  multi-zone,  upwind,  characteristic-based  code,  which  solves  the  integral 
form  of  the  Reynolds-Averaged  Navier-Stokes  equations.  GASP  models  the  fluid  continuum  with  a 
finite  number  of  control  volumes  and  incorporates  general  thermo-chemistry  modeling  with  an  ex¬ 
tensive  thermo-chemical  database.  The  latest  version,  GASP  v4,  has  more  efficient  time-integration 
schemes  including  implicit  treatment  of  zonal  boundaries.  GASP  v4  uses  a  distributed-memory  par¬ 
allel  architecture  in  order  to  exploit  a  wide  variety  of  platforms.  These  schemes  greatly  improve  the 
efficiency  of  the  code,  reducing  the  amount  of  time  needed  to  perform  a  computation.  AeroSoft’s 
SENSE  software  is  based  on  the  continuous  sensitivity  equation  method  and  provides  flow-field  sen¬ 
sitivities  to  a  wide  array  of  flow  and  geometric  design  variables.  SENSE  is  available  as  a  commercial 
package  and  has  been  used  in  a  variety  of  applications  of  aerospace  interest  [5,  6]. 
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The  initial  modeling  effort  included  modifications  to  GASP’s  thermo-chemical  database,  im¬ 
provements  to  the  species  diffusion  model,  and  modifications  to  account  for  laser  power  extraction. 
The  thermo-chemical  database  manager  has  been  modified  to  allow  chemical  species  to  be  present  at 
multiple  energy  states,  e.g I  (2Pi/2)  and  I  {^^3/2)  5  a  feature  characteristic  of  COIL  chemistry  mod¬ 
els.  A  20-reaction,  10-species  COIL  chemistry  model  has  also  been  added  to  the  thermo-chemical 
database.  The  species  diffusion  model  in  GASP  v4  has  been  extended  to  account  for  effective  dif¬ 
fusion  coefficients  as  well  as  the  pressure-driven  diffusion  effects  that  are  important  in  COIL  lasers. 
The  diffusion  models  have  been  tested  for  injection  problems  involving  helium  and  iodine  at  con¬ 
ditions  similar  to  the  RADICL  injector.  Software  modifications  have  also  been  made  to  GASP  v4 
to  account  for  laser  power  extraction  in  the  optical  cavity.  This  task  involved  adding  lasing  source 
terms  to  the  iodine  species  continuity  equations  and  the  global  energy  equation,  implementing  a 
geometric  optics  resonator  model,  and  including  conservative  interpolation  methods  between  the 
flow-solver  grid  and  the  optics  grid.  Power-On  results  have  been  obtained  for  the  RADICL  device. 

Numerical  studies  [10]  have  shown  that  water-vapor  condensation  in  the  COIL  configuration  may 
have  a  significant  impact  on  power  output.  A  homogeneous  water- vapor  condensation  model  similar 
to  that  described  in  Ref  [10]  has  been  developed  and  implemented  into  GASP  v4,  and  results  have 
been  obtained  for  the  RADICL  device.  Wall  catalysis  effects  also  reduce  COIL  efficiency,  as  excited 
species  contact  the  nozzle  wall  and  become  deactivated.  A  wall  catalysis  model  similar  to  that 
described  in  Ref  [8]  has  been  developed  and  implemented  into  GASP  v4.  The  model  includes  wall 
deactivation  and  recombination  reactions  and  is  consistent  with  both  the  10-specie  COIL  chemistry 
model  and  the  effective  diffusion  model.  Results  with  wall  catalysis  effects  have  been  obtained  for 
the  RADICL  system. 

Work  toward  the  second  goal,  i.e.,  (2)  above,  has  focused  on  extending  the  capability  of  SENSE  to 
include  models  important  for  COIL  systems.  Modeling  enhancements  to  SENSE  include  the  addition 
of  the  COIL  chemistry  model,  the  multicomponent  diffusion  model,  LeRC  thermodynamic  curve  fits, 
and  modules  to  compute  the  gain  and  gain  sensitivity.  The  capability  to  compute  flow  sensitivities 
to  Arrhenius  reaction  rates  has  also  been  implemented  into  SENSE.  A  sensitivity  analysis  module 
has  also  been  developed  for  the  laser  optics  resonator  and  is  based  on  differentiation  of  the  geometric 
ray-tracing  algorithm. 

Work  toward  the  third  goal,  i.e.,  (3)  above,  has  been  focused  on  developing  the  software  to 
couple  the  flow-field  sensitivity  module,  SENSE,  and  the  laser  optics  resonator  sensitivity  module. 
Numerical  studies  indicate  that  a  Gauss-Seidel  blocking  of  the  sensitivity  equations  permits  the  re¬ 
use  of  single-discipline  sensitivity  procedures  in  such  problems  with  multi-disciplinary  physics.  The 
COIL  sensitivity  package  has  been  tested  on  simple  geometries,  for  both  cold  flow  and  flow  with  the 
laser  on. 


2  Modifications  to  the  Thermo-Chemical  Database 

The  COIL  laser  operates  on  the  2(Pi/2)  2(^3/2)  transition  between  the  spin  orbit  components 

of  atomic  iodine.  This  transition  is  “pumped”  through  energy  transfer  from  the  excited  02(*A) 
state  of  molecular  oxygen.  Continuous  wave  operation  is  achieved  by  injecting  molecular  iodine  into 
a  primary  flow  containing  the  excited  states,  C^A)  and  02(:£),  and  the  ground  state,  02(3£), 
of  molecular  oxygen.  Once  the  molecular  iodine  has  been  injected  into  the  flow,  energy  pooling 
processes  involving  the  excited  states  of  oxygen  result  in  the  dissociation  of  molecular  iodine,  I2, 
into  atoms.  This  is  followed  by  rapid  energy  transfer  and  the  establishment  of  a  population  inversion 
in  atomic  iodine  according  to  the  “pumping”  reaction 

1  +  02(1A)^r  +02(3S),  (1) 

and  the  various  competing  mechanisms  that  work  to  reduce  I*.  Note  that  in  Eq  (1)  I  and  J*  denote 
the  2(P^/2)  and  2(Pi/2)  states  of  atomic  iodine  respectively.  During  the  lasing  process,  it  is  the 
excited  I*  state  of  iodine  that  is  stimulated  to  emit  photons,  while  the  energy  in  this  system  is 
stored  in  the  02 (2  A)  state.  Because  of  the  relatively  large  concentration  of  02 (x  A),  each  iodine 
atom  can  be  repumped  and  cycled  many  times  before  leaving  the  lasing  cavity. 
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No. 

Reaction 

Mechanism 

Kf(Kg  -  moZe,  m3,  s) 

1 

02(x A)  +  CM1  A)  -4  O2CS)  +  02(3S) 

I2  Dissociation  (+02(1£)) 

1.626  x  104  1 

2 

02(1S)  +  H2O  -4  02?A)  +  h2o 

I2  Dissociation  (—  02(1^)) 

4.035  x  109 

3 

02? A)  +  02(3s)  -4  02(3E)  +  02(3E) 

Competing  (-02(1A)) 

9.635  x  102 

4 

02?  A)  +  H20  -4  02(3£)  +  H20 

Competing  (—  O^1  A)) 

2.409  x  103 

5 

02?A)  +  C12  -4  02(3Z)  +  C12 

Competing  (-02(1A)) 

3.613  x  103 

6 

02?A)  +  He  -4  02(3S)  +  He 

Competing  (-02(1A)) 

4.818  x  10° 

7 

I2  "1"  O2  (*  2)  -+  2  7  +  O2  (3  2) 

I2  Dissociation  (+  7) 

2.409  x  109 

8 

/2  +  02  0  5j)  -+  I2  +  02(35j) 

72  Dissociation  (—  02  (XS)) 

9.635  x  109 

9 

I2  +  02  C  A)  — >  /I  +  02  (3S) 

72  Dissociation  (+72) 

4.215  x  106 

10 

h  + r  -4  /2*  + 1 

Competing  (—7*) 

2.288  x  1010 

11 

I2  +  02  (XA)  -4  2 1  +  02  (3  X) 

I2  Dissociation  (+7) 

1.807  x  1011 

12 

/|  +  02(3S)  -4  I2  +  02(3E) 

72  Dissociation  (— 72) 

3.011  x  1010 

13 

j2*  +  h2o  -4  /2  +  h2o 

I2  Dissociation  (— 7|) 

1.807  x  1011 

14 

I 2  4"  He  — y  I 2  4"  He 

I2  Dissociation  (-72) 

2.409  x  109 

15 

I  +  02?A)  ^7*  +  02(3S) 

Pumping  (+7*) 

4.697  x  1010 

16 

I  +  02?A)  -4  /  +  02(3s) 

Competing  (-02(1  A)) 

6.022  x  105 

17 

I* +  02? A)  -+I  +  O2? E) 

Competing  (-7*) 

6.022  x  107 

18 

I*  +  02?A)  -4  J  +  02(xA) 

Competing  (—7*) 

6.624  x  107 

19 

I*  +  I  -4  21 

Competing  (-7*) 

9.635  x  106 

20 

I*  +  H20  -4  /  +  H20 

Competing  (-7*) 

1.204  x  109 

Table  1:  10-specie,  20-reaction  COIL  chemistry  model 


Species 

Database  Name 

02(3S) 

0_2 

CM1  A) 

(0_2)_[sd] 

02?Z) 

(0_2)_[ss] 

h 

1.2 

n 

(I_2)_[star]  . 

I(2P3/  2) 

I 

I(2Pin) 

(I) _ [star] 

Table  2:  Thermo-chemical  database  species  naming  convention. 


During  this  Phase-I  effort  the  COIL  chemistry  model  in  Table  1  has  been  added  to  GASP’s 
thermo-chemical  database.  The  7 2  dissociation  mechanism  consists  of  the  two  separate  reaction 
paths  described  by  reactions  7  and  11.  Reactions  1,  2,  8,  9,  and  12  —  14  affect  the  I2  dissociation 
process  through  the  production  (or  destruction)  of  02  (x£)  and  7 2.  Reaction  15  is  the  pumping 
reaction,  and  reactions  3-6,  10,  and  16  -  20  work  to  deactivate  the  system  by  reducing  either 
O2OA)  or  7*.  The  forward  reaction  rates  for  each  reaction  are  constant  and  given  in  Table  1. 
The  backward  rates  are  computed  using  the  forward  rates  and  the  equilibrium  constant  which  is 
determine  through  the  LeRC  curve  fits  and  the  minimization  of  Gibb’s  free  energy. 

As  indicated  in  the  above  discussion,  the  COIL  chemistry  mechanism  is  used  (in  part)  to  model 
transitions  in  the  internal  energy  states  of  the  oxygen  and  iodine  species.  As  a  result,  the  various 
energy  states  of  a  given  molecule  or  atom  must  be  represented  as  a  separate  species  in  the  thermo¬ 
chemical  database.  To  accommodate  this  feature,  a  tagging  convention  has  been  added  to  the 
thermo-chemical  database  manager  which  allows  a  chemical  species  to  be  present  at  multiple  energy 
states.  The  various  energy  levels  of  oxygen  and  iodine  have  been  entered  into  the  database  as  shown 
in  Table  2.  Figures  1  and  2  illustrate  the  new  naming  convention  as  shown  in  the  species  and  model 
windows  of  the  database  manager  for  the  COIL  mechanism. 
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Figure  1:  Species  window  illustrating  new  species  naming  convention  and  species  parsing 


Figure  2:  Model  window  illustrating  the  new  COIL  chemistry  model 


3  Effective  Diffusion  Model 


Pressure  diffusion  terms  have  been  shown  to  be  important  in  the  injector  region  of  the  COIL  where 
a  2:1  or  greater  pressure  ratio  must  exist  for  sonic  injection  to  be  achieved.  During  Phase-I  an 
effective  diffusion  model  with  pressure-gradient  terms  has  been  added  to  the  viscous  formulation 
in  GASP  v4.  The  diffusion  terms  appear  in  the  species  continuity  equations  as  well  as  the  global 
energy  equation.  The  continuity  equation  for  species  s  is  given  by 


d_ 

dt 


ps  dV  + 


V  •  n)  dA  = 


)dA  +  JIJ 


uic  dV . 


(2) 


where  ps  is  the  density  of  species  s,  V  is  the  velocity  vector,  us  represents  the  chemistry  production 
term  for  species  s,  and  Vs  is  the  diffusion  velocity  vector  for  species  s.  Prior  to  this  contract  the 
species  diffusion  velocity  was  given  by  Fick’s  law 


Ps^s  —  pDsV(Ps/p)  ) 


(3) 


with  a  species  diffusion  coefficient  determined  using  a  constant  Schmidt  number  as 


D8  = 


P 

pSc 


(4) 


The  effective  diffusion  model  represents  a  modeling  improvement  in  that  it  accounts  for  individual 
species  diffusion  coefficients  as  well  as  pressure-driven  diffusion  effects.  The  effective  diffusion  model 
for  species  s  is  given  by 


PsVjs  —  Tt  Dsm 


+ 


(*-?) 


Vp 

p 


N 


+  7 

P  3= 1 


VXj  + 


(5) 


where  Ms  is  the  species  molecular  weight,  Dsm  is  the  species  effective  diffusion  coefficient,  7$  is  the 
mixture  concentration,  and  Xs  is  the  species  mole  fraction.  Definitions  for  the  species  concentration, 
7s,  the  mixture  concentration,  and  the  species  mole  fraction  are  given  respectively  as 


> 


N 

It  =  ^7^» 
$=1 


(6) 


The  species  effective  diffusion  coefficients  are  given  in  terms  of  a  mixture  summation  of  the  binary 
diffusion  coefficients,  Vij,  as  follows 


-D  0771  — 


1  -Xs 


Ef=ijVs  X:/Vsj 


(7) 


and  satisfy  the  mass  conservation  constraint  that 


N 

s=l 


(8) 


Closure  of  this  system  of  equations  requires  relations  for  the  binary  diffusion  coefficients,  T>ij ,  which 
have  been  determined  by  Hirschfelder  [7]  using  appropriate  collision  integrals.  The  resulting  equa¬ 
tions  are 


Vij  =  0.018828  T3/2 


pafjUij 


(9) 


where 


<Tii  — 


CFi  +  Oj 


(10) 


5 


'  0.75571  0.68191  ' 

y.*0.08742  ■*"  j.*0. 84910  ’ 


(11) 


T*  =  —  ,  (12) 

€ij 

and 

Cij  =  y/eiTj  .  (13) 

Actual  implementation  of  the  effective  diffusion  model  has  been  performed  by  separating  the 
effects  of  mole-fraction  gradients  and  pressure  gradients.  This  feature  allows  the  option  to  run  with 
or  without  the  pressure  diffusion  terms.  Separating  these  effects  in  Eq  (5)  yields 

N 

PsV}  =  -7t  Ms  Dsm  VXs  +7 Dim  VXi  ,  (14) 

P  3= 1 

and 


such  that 


P.V.  =  P.V*  +  P.V>  . 

Also  to  facilitate  ease  of  implementation,  the  following  variables  have  been  defined 

G  x  =  MsDsmVXs, 

and 

Gps  =  Ms  Ds 

Utilizing  these  definitions,  Eqs  (14)  and  (15)  become 


PsV}  =  -It 


(16) 

(17) 


(18) 


(19) 


and 


PtV*  =  —7 1 


N 


Vp 

P 


(20) 


The  diffusion  term  resulting  from  mole-fraction  gradients,  Eq  (19),  behaves  in  a  similar  fashion 
to  Fields  law  and  results  in  molecular  diffusion  away  from  mole-fraction  gradients.  The  pressure 
diffusion  term,  Eq  (20),  results  in  molecular  diffusion  of  “lighter”  species  away  from  pressure  gradi¬ 
ents  and  “heavier”  species  toward  pressure  gradients.  To  see  this  we  consider  a  mixture  of  helium 
and  iodine  with  molecular  weights  Mne  =  4.003  and  M/2  =  253.82  respectively.  The  species 
densities,  pressure  and  temperature  are  taken  as  pHe  —  0.010395  Kg/m3,  pi2  =  0.007605  Kg/mz, 
p  =  9576  N/m2  and  T  =  438.5  K.  These  conditions  results  in  mass  fractions  of  pHe/p  =  0.5775  and 
p/2/ p  ~  0.4225  and  mole  fractions  of  XHe  =  0.988564  and  of  Xh  —  0.011406.  Note  that  the  term 
(Xs  —  Ps/p)  in  Eq  (18)  will  be  positive  for  a  species  whose  mole  fraction  is  greater  than  its  mass 
fraction,  indicating  a  relatively  “light”  species.  This  same  term  will  be  negative  for  a  species  whose 
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“heavy”  species.  Evaluating  the 


(21) 


(22) 

direction  opposite  the  pressure- 
gradient  vector,  while  the  heavier  iodine  molecules  diffuse  in  the  same  direction  as  the  pressure- 
gradient  vector.  The  net  effect  is  for  light  and  heavy  particles  to  tend  to  separate  in  regions  of  large 
pressure  gradients. 

3.1  Helium-Iodine  Injection  into  an  Infinite  Chamber 

As  a  first  test  problem  for  the  new  diffusion  model,  we  considered  the  two-dimensional,  supersonic 
injection  of  helium  and  iodine  into  an  infinite  chamber.  Figure  3(a)  shows  the  81  x  161  computational 
mesh  used  for  this  case.  The  injector  face  extends  from  y  =  0.0  m  to  y  =  ±0.0005  m  with  flow 
entering  the  domain  from  left  to  right.  Solid  chamber  walls  extend  above  and  below  the  injector 
face.  Grid  points  are  clustered  in  the  horizontal  direction  near  the  wall  and  injector  face  and 
in  the  vertical  direction  near  the  injector  side  walls.  The  jet  Mach  number  was  M  =  1.3,  the 
temperature  was  T  —  288  AT,  and  the  pressure  was  p  =  28728  N/m2.  The  species  densities  were 
pHe  —  0.047480  Kg/mz  and  pj2  =  0.034737  Kg/m3.  Initial  conditions  in  the  chamber  consisted 
of  pure  helium  at  Mach  number  of  M  =  0.01,  a  temperature  of  T  =  300  AT,  and  a  pressure  of 
p  =  957 6  N/m2.  The  initial  species  density  for  helium  was  pne  —  0.015368  Kg/m3.  The  jet 
conditions  were  imposed  using  a  split-flux  free-stream  boundary  condition  and  the  chamber  wall 
was  treated  using  a  no-slip  wall  boundary  condition.  All  other  boundaries  were  treated  using  first- 
order  extrapolation  from  the  interior. 

Figures  3(b),  3(c),  and  3(d)  show  the  resulting  Mach  number  contours  partially  superimposed 
with  streamtraces,  pressure  contours,  and  iodine  mass  fraction  contours.  Figures  3(b)  and  3(c) 
indicate  that  the  helium-iodine  mixture  undergoes  a  strong  expansion  as  it  leaves  the  injector, 
resulting  in  the  lower  pressures  and  higher  Mach  numbers  shown  by  the  contours  plots.  The  very 
strong  pressure  gradients  near  the  injector  side  walls,  i.e.,  y  =  ±0.0005 m,  result  in  large  pressure- 
diffusion  effects  which  work  to  separate  the  helium  from  the  iodine.  This  can  be  seen  in  Fig  3(d), 
where  the  iodine  concentrations  in  the  near- wall  region  are  essentially  zero.  Comparison  of  the 
streamlines  in  Fig  3(b)  and  the  iodine  mass  fraction  contours  in  Fig  3(d)  indicate  that  the  diffusion 
process  occurs  in  a  direction  transverse  to  the  mean  flow  velocity.  Figure  4  shows  the  normalized 
pressure  profile  and  helium  and  iodine  mass-fraction  profiles  for  the  vertical,  i  =  10  grid  line.  This 
grid  line  is  very  close  in  proximity  to  the  injector  face  and  chamber  walls.  The  normalized  pressure 
profile  shows  a  large  pressure  gradient  at  y  =  +.0005  m  near  the  injector  side  wall.  The  iodine  mass 
fraction  profile  shows  a  maxima  in  the  high  pressure-gradient  region,  while  the  helium  mass  fraction 
profile  shows  a  minima.  These  results  are  consistent  with  the  pressure-diffusion  effects  previously 
described.  As  helium  molecules  encounter  the  pressure  gradient  they  tend  to  diffuse  in  the  +y 
direction  away  from  the  pressure  gradient,  resulting  in  a  minima.  Iodine  molecules  tend  to  diffuse 
in  the  —  y  direction  toward  the  pressure  gradient  resulting  in  the  observed  maxima. 

3.2  Transverse  Helium-Iodine  Injection 

As  a  second  test  problem  we  considered  the  two-dimensional,  transverse  injection  of  helium  and 
iodine  into  a  pure  helium  free  stream.  Figure  5(a)  shows  the  113  x  81  computational  mesh  used 
for  this  case.  The  lower  surface  is  composed  of  a  solid  wall  with  the  injector  face  located  between 
x  =  0.0075  m  and  x  =  0.0077  m.  Grid  points  are  clustered  in  the  vertical  direction  near  the  wall 
and  injector  face  and  in  the  horizontal  direction  near  the  injector  side  walls.  The  jet  conditions 


mole  fraction  is  less  than  its  mass  fraction,  indicating  a  relatively 
terms  in  Eq  (20)  for  the  specified  conditions  yields 

PHeVpHe  =  -1.2456  x  10-4  ^  , 

p 

and 

Phvph  =  +1.2456  x  lO"4  ^  . 

These  results  show  that  the  lighter  helium  molecules  diffuse  in  a 
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(a)  Two-Dimensional,  81  x 
161  grid. 
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(c)  Pressure  contours. 


(d)  Iodine  mass  fractions. 


Figure  3:  Helium-Iodine  injection. 
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Helium  and  Iodine  Mass  Fractions 
M=1.3,  He-I2  Injection  Case  (1=10) 


Figure  4:  Helium  and  iodine  mass  fractions  (station  1=10) 


are  identical  to  those  in  Section  3.1.  The  free-stream  conditions  are  identical  to  the  initial  chamber 
conditions  in  Section  3.1,  except  the  Mach  number  is  M  =  0.8  with  the  flow  direction  being  left  to 
right.  As  before,  the  jet  conditions  were  imposed  using  a  split-flux  free-stream  boundary  condition 
and  the  lower  wall  was  treated  using  a  no-slip  wall  boundary  condition.  The  subsonic  inflow  boundary 
was  treated  by  fixing  the  all  the  flow  variables  except  pressure  to  the  free-stream  values.  The  pressure 
at  the  inflow  boundary  was  extrapolated  from  the  interior.  All  other  boundaries  were  treated  using 
first-order  extrapolation  from  the  interior. 

Figure  5(b)  shows  the  iodine  mass  fraction  contours  for  the  complete  domain.  Figure  6(a)  shows 
a  close-up  view  of  the  pressure  contours  and  streamtraces  in  the  injector  region.  Comparison  of  these 
figures  shows  three  distinct  flow  regions;  the  primary  flow  of  helium  which  expands  around  the  jet 
to  supersonic  speeds,  the  jet  flow  which  contains  regions  of  iodine  concentrations  that  are  slightly 
higher  than  jet  free-stream  values,  and  a  recirculation  region  just  to  the  right  of  the  injector  with 
lower  than  jet  free-stream  concentrations  of  iodine.  The  variation  in  iodine  concentrations  between 
the  jet  flow  region  and  the  recirculation  region  is  caused  by  the  pressure  diffusion  effects  that  tend 
to  separate  iodine  and  helium.  Figure  6(b)  shows  the  normalized  pressure  profile  and  iodine  mass 
fraction  profiles  for  the  j  =  2  grid  line  in  the  injector  region.  The  normalized  pressure  profile 
shows  two  large  pressure-gradient  regions  in  the  vicinity  of  the  injector  side  walls.  The  iodine  mass 
fraction  profile  (for  the  effective  diffusion  model)  shows  peaks  at  these  locations,  indicating  that  the 
iodine  molecules  are  diffusing  in  the  direction  of  the  pressure-gradient.  This  causes  the  jet  flow  to 
have  slightly  higher  concentrations  of  iodine  than  expected.  At  the  same  time  helium  molecules  are 
diffusing  away  from  the  pressure  gradient  causing  the  recirculation  region  to  be  composed  of  primarily 
helium.  Figure  6(b)  also  shows  the  iodine  mass  fraction  profile  for  a  simulation  using  Fick’s  law, 
which  has  no  pressure  terms.  In  this  case  both  the  jet  and  recirculation  region  have  helium  and 
iodine  concentrations  consistent  with  the  jet  free-stream  conditions.  Figures  7(a)  and  7(b)  show 
iodine  mass  fraction  contours  in  the  vicinity  of  the  injector  for  both  the  effective  diffusion  model 
and  Fick’s  law. 

4  Coupling  GASP  with  a  Laser  Optics  Resonator  Model 

COIL  analysis  codes  currently  in  use  at  AFRL/DE  can  be  viewed  as  coupled  fluid  dynamics  and 
Laser  Optics  Resonator  (LOR)  codes.  The  MINT  implementation  [1]  of  the  COIL  simulation  couples 
a  pde-based  model  of  the  mechanics  for  compressible,  viscous,  chemically  reacting  flow  with  a  ray- 
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(a)  2-Dimensional,  113  x  81  grid. 


(b)  Iodine  mass  fraction  contours. 

Figure  5:  Helium-iodine  transverse  injection  case. 

tracing  model  [2]  for  the  optics  resonator.  The  GASP  implementation  described  in  this  section 
follows  the  same  construct. 

4,1  Lasing  Source  Terms 

The  CFD/LOR  coupling  occurs  through  laser  power  extraction  terms.  Power  extraction  in  the  lasing 
cavity  manifests  itself  in  the  form  of  source  terms  in  the  iodine  species  continuity  equations,  and 
in  the  global  energy  equation.  The  species  continuity  equations  are  affected  because  components  at 
different  energy  states  are  represented  as  separate  species.  In  a  COIL  laser,  iodine  in  the  excited, 
I(2Pi/2),  energy  state  is  stimulated  to  emit  photons  (i.e.,  stimulated  emission).  Upon  emission  of 
the  photon,  the  iodine  atom  assumes  the  lower,  I  (2P3/2) ,  energy  state.  As  a  result,  power  extraction 
affects  the  concentration  of  these  two  states.  Power  extraction  also  accounts  for  a  net  energy  loss 
(through  emission  of  photons)  to  the  flow.  The  modified  species  continuity  equation  is  given  as 

Ft  III  P'dV  +  ii)  dA  =  fA(-p8V,  n)  dA  +  jjj  uadV±  jjj  ^jj^dV,  (23) 

where  the  last  term  on  the  right-hand  side  of  Eq  (23)  is  the  power-extraction  source  term.  The 
source  term  is  composed  of  the  gain,  a,  the  two-way  average  intensity,  /,  the  molecular  weight,  Af#, 
and  the  energy  of  a  photon,  hu  =  90956  J/g  -  mole .  The  field  variable,  a,  represents  the  gain  and 
is  defined  as 

(24) 

where 
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(a)  Pressure  contours  and  streamtraces. 


Iodine  Mass  Fractions 
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(b)  Normalized  pressure  and  iodine  mass  fractions. 


Figure  6:  Helium-iodine  transverse  injection  case  (close-up  of  injector). 
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(b)  Iodine  mass  fractions  for  Fick’s  law. 

Figure  7:  Helium-iodine  transverse  injection  case  (close-up  of  injector). 
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Here  TV/*  and  IV/  represent  number  densities  for  the  excited  and  ground  states  of  atomic  iodine 
and  <f>(v)  is  the  Voight  line-shape  function.  The  constants  A  and  A  are  taken  as  A  =  5.0  sec~l  and 
A  =  1.31527  x  10“4  cm.  The  plus  sign  for  the  lasing  source  term  in  Eq  (23)  is  used  for  the  ground 
state  of  iodine  and  the  minus  sign  is  used  for  the  excited  state.  No  other  species  in  the  model  has 
a  source  term  contribution  from  power  extraction.  The  global  energy  equation  has  a  similar  source 
term,  —a  /,  to  account  for  the  energy  removed  from  the  flow. 

The  two-way  average  intensity,  J,  is  provided  as  output  from  the  LOR  model  (in  this  work,  the 
ray-trace  model  of  Ref  [2]),  and  is  a  function  of  the  gain,  the  mirror  radii  and  reflectivities,  the 
geometry  of  the  optical  cavity,  and  various  transmission  and  diffractive  losses. 


4.2  The  Coupled  CFD/LOR  Model 

In  the  coupled  CFD/LOR  model  we  have  two  systems,  a  fluid-dynamic  system  denoted  by  Tlf  and 
a  laser  optics  resonator  system  denoted  by  7 lp 

Kf(ufiup)  =  0  (29) 

7 lp{iif,up)  =  0.  (30) 

Here,  the  variable  u/  describes  the  fluid-dynamic  state  and  the  variable  up  describes  laser-power 
state.  An  overview  of  the  coupled  CFD/LOR  analysis  procedure  is  as  follows: 

1.  The  CFD  flow  solver  uses  a  pseudo-time  iteration  to  drive  the  flow  residual  ( i.e .  Tlf)  to  zero. 
The  state- variable  from  the  laser-power  extraction  model,  generically  up,  represents  the  two- 
way  average  intensity,  I.  The  intensity  appears  in  the  fluid-dynamic  system  in  the  continuity 
source  terms  for  the  iodine  species  and  in  the  global  energy  equation  as  described  above. 

2.  For  a  fixed  intensity  field  I  the  flow  iteration  proceeds  to  reduce  the  flow-residual.  The  flow- 
state  variable,  generically  u/,  includes  the  flow  velocities,  thermodynamic  variables  and  species 
densities.  From  these  known  quantities  we  compute  a  gain  field,  a,  using  Eq  (24).  Values  in 
the  gain  field  are  computed  pointwise  from  the  flow-field  values. 

3.  The  gain  field  is  passed  to  the  ray-tracing  module  for  the  calculation  of  an  updated  intensity 
field.  There  are  several  steps  in  this  linking  process  that  arise  from  modeling  choices.  In 
particular,  reduced  flow  simulations  are  typically  based  on  a  unit-cell  approximation  that 
imposes  a  periodic  variation  of  the  flow- variables  in  the  direction  of  the  optical-axis.  As  a 
result,  the  gain  field  is  averaged  in  the  direction  of  the  optical  axis  across  the  active  region. 
The  average  gain  field  from  the  flow-solver  mesh  is  then  interpolated  to  a  two-dimensional 
polar  mesh  consistent  with  the  ray-trace  algorithm.  This  interpolated  gain  is  then  passed  to 
the  ray-trace  algorithm. 

4.  Steady-state  lasing  occurs  for  the  ray-trace  algorithm  when  the  round-trip  gain-equals-loss 
condition  is  satisfied.  Such  a  steady  field  is  achieved  only  when  gain  production  and  laser 
extraction  are  in  equilibrium.  As  a  result,  the  flow  and  power  models  must  be  solved  simulta¬ 
neously.  In  order  to  make  sense  of  a  sequential  iteration  procedure  one  must  define  a  unique 
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Figure  8:  Schematic  of  the  geometric  optical  model 


intensity  field  that  is  associated  with  the  given  gain  field.  The  power  extraction  module  based 
on  [2]  does  this  by  calculating  a  return  intensity  and  using  an  under-relaxation  update.  This 
process  is  described  in  more  detail  in  the  following  section. 

5.  The  updated  intensity  field  from  the  power  extraction  module  is  communicated  back  to  the  flow 
model.  These  intensities  are  then  interpolated  from  the  polar  optics  mesh  to  the  flow-solver 
mesh  where  they  are  used  to  compute  updated  lasing  source  terms. 

4.3  Geometric  Optics  Ray- Tracing  Model 

The  main  work  of  the  laser  optics  resonator  model  is  done  within  the  ray-tracing  algorithm  [2].  A 
significant  code  fragment  is  reproduced  here. 

1000  continue 
c 

c  calculation  of  intensity  on  subsequent  iterations,  i.e.,  iopt  =  0 
c 

sdifmax  =0.0 
do  1500  i  =  l,nthetal 
do  1200  j  =  l,nrmax 

rod(j,i)  =  RABS(rmaxls(i)  *  (xdist  -  xa(j))  /  distm) 
sum  =  storg(j,i)  -  scripl 
alphaq(l)  =  storg(j,i) 
do  1100  iq  =  2,m 

alam  =  SQRT(1.0  +  (rod(j,i)  *  afl(iq))**2) 
alphaq(iq)  =  xratio(j,iq)  *  (storg(iset (j , iq)-l, i)  - 
1  storg(iset ( j , iq)  , i) )  +  storg(iset (j , iq)  , i) 

sum  =  sum  +  alam  *  (alphaq(iq)  -  scripl) 

1100  continue 

s(j)  =  2.0  *  glength  *  sum 
exps  =  EXP(RMIN(expmin,s(j))) 

c  ensure  that  stori(j,i)  is  not  changed  on  restart 

if(ifirst  .eq.  1)  then 
iret(j)  =  inot(j,i) 
else 

iret(j)  =  inot(j,i)  *  rlr2m  *  exps  *  prodj(j) 
endif 

inot(j,i)  =  coninot  *  iret(j)  +  (1.0  -  coninot)  *  inot(j,i) 
expl  =  EXP (RMIN(expmin, glength  *  (alphaq(l)  -  scripl))) 
terml  =  (expl  -  1.0)  /  (glength  *  (alphaq(l)  -  scripl)) 
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c 


1200 

c 

1500 


term2  =  1.0  +  prodj(j)  *  (rlr2m  *  exps  /  (rml  *  expl))  / 
(1.0  -  del(j)) 

stori(j,i)  =  inot(j,i)  *  terml  *  term2 

(several  lines  omitted) 
cont inue 

(several  lines  omitted) 
continue 


The  i  and  j  loops  are  stepped  through  the  gain  region  using  polar  coordinates  to  span  the 
(x,y)  -  plane  in  a  coordinate  system  that  is  convenient  for  the  geometry  of  the  optical  cavity.  Each 
(i ,  j)  pair  labels  a  point  on  the  out-coupling  mirror  and  the  calculations  within  the  loop  implement 
Eq  (7)  in  [2,  page  5].  Figure  8  shows  a  top- view  schematic  of  the  optical  cavity.  Each  grid  point 
on  the  polar  optics  mesh  labels  a  point  on  the  out-coupling  mirror.  A  ray  of  initial  intensity,  Jo, 
leaving  a  point  on  the  flat  out-coupling  mirror  in  an  orthogonal  direction  will  return  (orthogonally) 
to  the  same  point  after  2m  passages  through  the  gain  medium.  After  each  passage  there  is  a  specular 
reflection  from  the  opposite  mirror.  The  return  intensity,  Ir,  is  given  by: 


Ir(J)  =  Io(J)(RiR2)m  es 


m 


m 


JJ  (1  Sp)s  —  2 Lg  Yjap 
p= i  p= i 


-£)  A, 


(31) 


where  Rx  and  R2  are  the  reflectivities  for  the  flat  outcoupler  and  spherical  mirror,  and  ap  is  the 
average  gain  along  the  p- th  ray  segment.  The  initial  intensity  field,  Jo,  for  the  next  invocation  of 
this  routine  is  computed  in  the  line: 

Jo  =  cu  Ir  H-  (1  —  oj)  Jo,  (32) 


or 


inot(j,i)  =  coninot  *  iret(j)  +  (1.0  -  coninot)  *  inot(j,i) 

Note  that  this  code  fragment  is  embedded  in  an  i-loop  so  the  i  variable  does  not  appear  explicitly  in 
iret  ( ) .  In  the  present  version  we  have  u  —  0.5  (coninot  =  0 . 5),  so  that  the  scheme  under-relaxes 
the  intensity  field.  When  the  coupled  system  converges  the  initial  and  return  intensity  values  will 
agree,  so  that,  in  effect,  the  relaxation  calculation  is  identical  to  Jo  =  Ir- 
Finally,  the  average  two-way'  intensity  is  computed  from 


t  _  /q(J)  [exp  [(a(J)  -  C)Lg]  -  1] 

HJ)-QW 


'  exp  [-  (a(J)  -  £)  Lg] 


(33) 


and  passed  to  the  CFD  flow  solver. 


4.4  Software  Modifications 

AeroSoft  has  performed  modifications  to  GASP  v4  to  facilitate  the  power-extraction  model.  The 
new  additions  to  GASP  include: 

1.  Memory  allocation  for  the  gain  and  intensity  fields. 

2.  Gain  computation  routine. 

3.  Lasing  source  term  and  Jacobian  routines. 

4.  Parallel  logic  to  gather  the  gain  from  zones  (within  the  lasing  cavity)  on  distributed  processors. 

5.  Parallel  logic  to  scatter  the  intensity  to  zones  (within  the  lasing  cavity)  on  distributed  proces¬ 
sors. 
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(a)  Bi-linear  interpolation  scheme. 


Figure  9:  Schematic  of  mapping  between  flow-solver  and  polar  optics  mesh. 


6.  Routines  to  define  the  optics  parameters  and  generate  the  polar  mesh. 

7.  Implementation  of  the  stable .  F  routine  for  the  ray-tracing  algorithm. 

8.  Interpolation  routines  between  the  flow-solver  mesh  and  the  polar  optics  mesh. 

9.  Iteration  strategy  for  updating  intensity. 

10.  Flow  visualization  capability  for  gain. 

5  Fully-Coupled  COIL  Simulations 

Initial  attempts  to  perform  a  fully-coupled,  power-on,  RADICL  simulation  yielded  solutions  with 
divergent  lasing  power.  The  resulting  flow  fields  exhibited  localized  regions  within  the  lasing  cavity 
where  the  intensity  continued  to  grow  unbounded.  We  hypothesized  that  this  instability  resulted 
from  an  “uncoupling”  of  grid  points  caused  by  the  bi-linear  interpolation  method  used  to  transfer 
the  gain  and  intensity  between  the  flow-solver  and  polar  optics  mesh.  Figure  9(a)  depicts  the  super¬ 
imposed  flow-solver  and  polar  optics  meshes.  Notice  for  the  particular  case  shown,  the  flow-solver 
mesh  is  very  coarse  relative  to  the  optics  mesh.  As  an  example,  we  consider  the  step  of  transferring 
the  intensity  from  the  optics  mesh  points  to  the  flow-solver  cell  center  (denoted  by  the  solid  black 
circle).  To  perform  this  step,  the  bi-linear  scheme  finds  the  four  closest  surrounding  cells  in  the  optics 
mesh  (t.e.,  the  squares).  These  four  points  are  then  used  in  a  bi-linear  interpolation  to  determine 
the  cell-center  value  for  the  intensity.  For  the  case  shown  in  Figure  9(a)  there  are  additional  mesh 
points  (the  triangles)  inside  the  flow-solver  cell  that  are  never  used  in  the  interpolation  procedure. 
The  intensity  at  these  “uncoupled”  points  has  no  influence  on  the  flow-field.  Points  similar  to  these, 
in  our  initial  RADICL  simulations,  coincided  with  locations  where  the  intensity  was  allowed  to  grow 
unchecked. 

In  order  to  eliminate  this  problem,  a  new  conservative  interpolation  scheme  has  been  developed 
to  transfer  the  gain  and  intensity  between  the  flow  solver  and  optics  meshes.  In  this  scheme,  a 
pseudo-cell  is  constructed  around  each  optics  grid  point.  The  intensity  at  a  flow-solver  cell  center 
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(34) 


would  then  be  influenced  by  each  of  the  overlapping  optics  cells  as  follows 

i  _  Ijk-Aji 
i_  ZjAAji  ■ 

Here,  the  index  i  represents  cells  in  the  flow-solver  mesh  and  the  index  j  represents  cells  in  the  optics 
mesh.  The  area  term  Aji  represents  the  area  of  overlap  for  cell  j  (of  the  optics  mesh)  onto  cell  i 
(of  the  flow-solver  mesh).  The  greyed  cell  in  Figure  9(b)  depicts  the  area  of  overlap  for  one  of  the 
optics  cells  onto  a  flow-solver  cell.  Using  this  method,  all  optics  mesh  points  would  be  included  in 
the  interpolation  stencil. 

5.1  Simple  Test  Case 

Power-on  simulations  have  been  performed  for  a  simplified  flow  geometry,  but  with  flow  conditions 
and  optical  parameters  consistent  with  the  RADICL  laser.  This  test  problem  reduced  the  debugging 
effort  and  served  to  demonstrate  the  coupled  processes  involved  in  the  new  power-extraction  model. 
Figure  10(a)  shows  the  grid  topology  for  the  151  x  51  flow-solver  mesh  and  the  101  x  51  polar 
optics  mesh.  The  number  of  grid  points  shown  in  this  figure  has  been  reduced  to  facilitate  viewing. 
Notice  that  the  x-direction  stretching  at  the  leading  and  trailing  edges  of  the  optics  mesh  is  nearly 
identical  to  the  flow  solver  mesh.  The  start  of  the  aperture  is  at  x  =  .06  m  and  the  aperture 
length  is  La  =  .05  m.  The  aperture  height  is  Ha  =  ,028  m.  The  distance  between  the  curved 
mirror  and  the  flat  out-coupling  mirror  is  Dm  =  0.8  m  and  the  radius  of  curvature  of  the  curved 
mirror  is  r 2  =  10.0  m.  The  coefficients  of  reflectivity  for  the  mirrors  are  Ri  =  0.90  and  R2  =  0.99 
for  the  out-coupler  and  curved  mirrors  respectively.  The  gain  length  was  taken  as  Lg  =  0.25  m. 
Distributed  losses  in  the  gain  medium  and  diffractive  losses  at  the  aperture  were  neglected.  The 
inflow  boundary  was  split  into  two  portions;  with  the  lower  stream  having  a  high  concentration  of 
I*  (high-gain  region),  and  an  upper  stream  with  very  little  I*  and  C^A)  (low-gain  region).  The 
low-gain  region  was  designed  to  account  for  the  region  of  low  gain  observed  near  the  upper  wall  in 
the  power-off  RADICL  simulations.  The  inflow  conditions  for  the  test  case  were  chosen  from  actual 
conditions  (upstream  of  the  lasing  cavity)  for  the  power-off  RADICL  simulation,  and  are  listed  in 
Table  3. 

For  this  test  case  the  uncoupling  problem  did  not  occur  and  the  power-on  solution  produced  a 
converged  value  for  the  lasing  power.  The  solid  line  in  Figure  11  shows  the  convergence  history 
for  the  lasing  power  in  this  case.  Figures  10(b)  and  10(c)  show  the  gain  and  intensity  contours 
for  this  simulation.  Values  for  the  inflow  gain  in  the  high  and  low-gain  regions  are  approximately 
a  =  1.065  1/m  and  a  =  0.00751  1/m  respectively.  The  maximum  gain  occurs  in  the  high-gain 
region  just  before  the  start  of  the  aperture  and  has  a  value  of  Oimax  =  1.308  1/m.  Past  the  start 
of  the  aperture,  the  gain  is  dramatically  reduced  through  interaction  with  the  lasing  source  terms. 
The  maximum  intensity  occurs  at  the  aperture  leading  and  trailing  edges  with  a  value  of  I  max  = 
3.6  x  109  W/rn2.  The  power  output  for  this  case  converged  after  approximately  1000  iterations  and 
was  P  =  13.1  KW. 

In  a  second  test,  the  same  case  was  run  but  with  a  coarsened  flow-solver  mesh.  The  coarse 
flow-solver  mesh  was  obtained  by  sequencing  five  levels  in  both  the  x  and  y  directions,  yielding  a 
31  x  11  grid.  The  101  x  51  optics  mesh  from  the  previous  test  above  was  used.  The  grid  mismatch 
in  this  case  should  increase  the  possibility  that  grid-point  uncoupling  similar  that  shown  in  Fig  9 
will  occur.  As  expected,  the  lasing  power  for  this  case  did  not  converge.  The  dashed  line  in  Fig  11 
shows  the  lasing  power  history  for  this  case. 

5.2  RADICL  Simulation 

Using  the  new  conservative  interpolation  scheme,  a  power-on  simulation  has  been  performed  for  the 
RADICL  device.  The  physical  domain  of  the  RADICL  simulation  is  depicted  in  the  coarsened  view 
of  the  computational  mesh  shown  in  Figure  13.  Figure  14(a)  shows  a  close-up  view  of  injector  region. 
Our  simulation  considered  only  the  upper  half  of  this  domain,  and  employed  symmetry  conditions 
about  the  centerline  of  the  nozzle.  In  the  lateral  direction  (z.e.,  into  the  page),  the  RADICL  device 
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(a)  Coarsened  flow-solver  mesh  and  polar  optics  mesh. 


(b)  Gain  contours  (1/m). 


(c)  Intensity  contours  (' Watts/m 2). 


Figure  10:  Test  laser  configuration  and  results. 
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Variable 

High-Gain  Region 

Low-Gain  Region 

PHe 

0.00260088  Wgjrn? 

0.000860848  Wgjm? 

Po2 

0.00208539  Kg/m3 

0.000811512  Kg/m3 

Po2p-  a) 

0.00161352  Kg/m3 

0.000831930  Kg/m3 

Po2(3e) 

2.29196  x  10~6  Kg/m3 

8.38014  x  10"8  Kg/m3 

Ph2o 

0.000247298  Kg/m3 

0.000110298  Kg/m3 

Ph 

0.000235656  Kg/m3 

7.09078  x  10~5  Kg/m3 

Pi 3 

8.97101  x  10-7  Kg/m3 

8.90463  x  10-9  Kg/m3 

Pi 

4.74918  x  10-5  Kg/m3 

1.79366  x  10"6  Kg/m3 

Pi* 

0.000319047  Kg/m3 

4.33212  x  10~6  Kg/m3 

Pci2 

0.00122394  Kg/m3 

0.000545616  Kg/m3 

u 

881 .8m/s 

990.9  m/s 

V 

0.0  m/s 

0.0  m/s 

w 

0.0  m/s 

0.0  m/s 

P 

1124.6  N/m2 

1200.0  N/m2 

Table  3:  Inflow  conditions  for  laser  test  case. 


Figure  11:  Lasing  power  vs.  iteration  number. 


consists  of  115  sets  of  injector  orifices  (on  both  the  upper  and  lower  walls),  each  set  consisting 
of  one  large  and  two  small  injector  orifices.  Figure  12  shows  three  such  sets  of  injectors.  Each 
injector  set  is  symmetric  about  the  centerline  plane  of  the  large  injector  orifice.  The  large  orifices 
have  a  diameter  of  Dl  =  0.08128  cm  and  the  small  orifices  have  a  diameter  of  Ds  =  0.04064  cm. 
We  now  make  a  key  assumption  that  the  bulk  effect  of  all  the  injector  sets  and  nozzle  side  walls 
has  a  minimal  influence  on  a  single  chosen  injector  set.  With  this  assumption  a  single  injector  set 
can  be  isolated  for  simulation  and  the  lateral  extent  of  the  domain  can  be  reduced  by  employing 
symmetry.  Our  simulation  was  performed  on  the  representative  sub-domain  of  the  COIL  flow  field 
termed  the  “unit  cell”  and  depicted  by  the  dashed  line  in  Figure  12.  The  unit  cell  consists  of  half 
of  a  large  orifice  and  one  small  orifice.  While  the  ramifications  of  this  approximation  are  not  fully 
understood,  our  goal  was  to  model  the  flow  field  with  a  sufficient  level  of  detail  so  as  to  capture  the 
physical  processes  necessary  to  test  our  coupling  procedures.  This  would  have  been  impossible  for 
the  complete  physical  domain  because  a  prohibitive  number  of  grid  points  would  be  required. 

The  RADICL  grid  system  consisted  of  four  zones;  a  6  x  21  x  9  mesh  for  the  large  injector,  a 
5  x  21  x  9  mesh  for  the  small  injector,  a  147  x  77  x  25  mesh  that  includes  the  upstream  portion  of 
the  mixing  nozzle  and  throat,  and  a  145  x  77  x  25  mesh  for  the  downstream  portion  of  the  nozzle 
that  contains  the  optical  cavity.  Grid  points  were  clustered  in  the  ^-direction  near  the  leading  and 
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trailing  edges  of  the  aperture  as  depicted  in  Figure  13.  A  total  of  564, 179  grid  points  were  used. 

Nozzle  inflow  conditions  are  specified  using  a  subsonic  inflow  boundary  condition  with  fixed  total 
pressure  and  temperature.  The  total  conditions  were  consistent  with  the  measured  static  values  for 
pressure  and  temperature  of  p  =  9961.7  N/m2  and  T  =  315  AT.  The  chemical  composition  at  the 
inflow  was  consistent  with  measured  molar  flow  rates  and  an  C^A)  yield  of  Y  =  0.41.  The  species 
inflow  densities  were  specified  as 


PHe  =  0.01117200  Kg/m3 
Po2(3e)  =  0.01544990  Kg/m3 
Po2(i  A)  =  0.01073640  Kg/m 3 
Po2{1'Z)  =  0.00000019  Kg/m3 
Ph2o  =  0.00139578  Kg/m3 
pci2  =  0.00828914  Kg/m3 . 


The  inflow  velocity  was  V  =  84.028  m/s.  The  injector  inflow  conditions  were  specified  using  the 
same  subsonic  inflow  boundary  condition.  The  total  temperature  was  taken  to  be  Tt  —  400  if,  and 
the  species  densities  were  specified  to  match  measured  molar  flow  rates  of  He  and  I %  and  are  given 
as  follows 


PHe  =  0.0246814  Kg/m3 
ph  =  0.0132298  Kg/m 3  . 

The  nozzle  wall  is  treated  as  a  T  =  315  if,  no-slip  wall. 

The  aperture  (and  associated  mirror  system)  was  located  between  x  =  0.08578  m  and  x  = 
0.14578  m.  The  length  of  the  aperture  was  La  =  0.06  m.  The  distance  between  the  spherical  mirror 
and  the  flat  out-coupling  mirror  is  Dm  =  0.8  m  and  the  radius  of  curvature  of  the  spherical  mirror  is 
r2  =  10.0  m.  The  coefficients  of  reflectivity  are  i?i  =  0.927  and  R2  =  0.995  for  the  flat  out-coupler 
and  spherical  mirrors  respectively.  The  gain  length  was  taken  as  Lg  =  0.254  m  (the  lateral  width 
of  the  nozzle).  Distributed  losses  in  the  gain  medium  and  diffractive  losses  at  the  aperture  were 
neglected.  The  polar  mesh  size  for  the  ray-tracing  algorithm  was  101  x  41  and  laser  optics  resonator 
model  was  called  on  every  third  flow-solver  iteration. 

The  solution  procedure  utilized  four  mesh  sequences  and  required  approximately  50  hours  on  a 
32  processor  SGI  Origin.  The  fine-grid  solution  required  1000  iterations  to  reduce  the  l2-n orm  of  the 
residual  approximately  2.5  orders  of  magnitude.  GASP  v4  required  10  GB  of  memory  to  execute  on 
the  fine  mesh  utilizing  fully-implicit  zonal  boundaries. 

Figure  14(b)  shows  the  Mach  number  contours  in  the  vicinity  of  the  injector  region  for  the  k  —  9 
nodal  plane.  The  Mach  number  contours  provide  some  indication  of  the  jet  penetration  into  the 
primary  flow.  The  flow  from  the  large  injector  reaches  the  nozzle  centerline  just  upstream  of  the 
throat.  Figure  15(a)  shows  the  averaged  gain  contours  for  the  downstream  section  of  the  nozzle 
including  the  optical  cavity.  Gain  values  are  the  highest  just  before  the  start  of  the  aperture,  with 
a  maximum  of  amax  =  1.062 1/m.  Past  the  start  of  the  aperture,  the  gain  is  dramatically  reduced 
through  interaction  with  the  lasing  source  terms.  The  experimentally  measured  gain  1  cm  upstream 
of  the  aperture  leading  edge  and  on  the  nozzle  centerline  was  a  =  0.80 1/m.  The  COIL  simulation 
predicted  a  gain  at  this  location  of  a  =  0.79941/m.  Figure  15(b)  shows  the  intensity  contours 
predicted  by  the  ray-tracing  algorithm  in  the  lasing  cavity.  The  maximum  value  of  intensity  occurs 
at  the  aperture  leading  and  trailing  edges  with  a  value  of  Imax  =  3.6  x  109  W/m 2.  The  laser  power 
output  predicted  for  this  simulation  was  P  =  7.524  KW .  The  experimentally  measured  value  of 
power  was  P  —  5.8  KW.  This  is  reasonable  agreement,  since  no  loss  terms  were  included.  Similar 
studies  ( e.g Buggeln  [1])  have  used  non-zero  values  for  mirror  scattering  and  diffractive  losses  to 
give  corrected  output  powers  that  are  up  to  20%  lower  than  the  zero-loss  power.  Assuming  similar 
behavior  for  this  case,  our  power  results  are  in  very  good  agreement  with  the  measured  value.  The 
results  of  this  study  were  presented  in  Eppard  [4]. 
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6  Water- Vapor  Condensation  Model 

Water  vapor  is  produced  as  a  by-product  in  the  COIL’s  singlet-delta  oxygen  generator.  Since  water 
vapor  is  very  efficient  at  deactivating  the  excited  state  of  iodine,  it  is  removed  before  the  singlet 
oxygen  is  mixed  with  iodine.  Unfortunately,  a  small  portion  of  water  vapor  escapes  the  removal 
system,  mixes  with  the  iodine  and  expands  through  the  nozzle.  The  rapid  decrease  in  temperature 
in  the  diverging  section  of  the  nozzle  results  in  the  nucleation  and  subsequent  growth  of  water 
droplets.  The  positive  aspect  of  condensation  is  that  it  reduces  the  amount  of  water  vapor  in  the 
mixture,  and  therefore  suppresses  the  deactivation  of  excited  iodine.  Unfortunately,  as  the  water 
vapor  condenses  the  latent  heat  release  into  the  flow  raises  the  temperature  considerably  in  the 
lasing  cavity.  This  increase  in  temperature  causes  a  significant  reduction  in  the  gain  and  therefore 
power  output. 

To  model  this  effect,  we  have  implemented  into  GASP  v4  a  homogeneous  water-vapor  conden¬ 
sation  model  similar  to  that  described  by  Masuda  [10].  The  model  includes  the  latent  heat  in  the 
global  energy  equation  and  accounts  for  inter-phase  mass  transfer  through  source  terms  in  the  water 
vapor  and  droplet  continuity  equations.  The  model  assumes  that  the  water  droplets  have  the  same 
velocity  as  the  gas  phase,  and  does  not  require  particle  tracking.  This  assumption  is  reasonable  since 
the  droplet  size  within  the  COIL  is  very  small.  Finally,  we  assume  that  the  water  droplets  are  in 
thermal  equilibrium  with  the  carrier  gas,  z.e.,  the  droplet  temperature  is  identical  to  the  gas-phase 
temperature.  These  assumptions  are  consistent  with  the  previous  COIL  analysis  of  Masuda  [10]. 

The  condensation  model  considers  the  liquid  phase  using  a  number  of  droplet  classes,  or  par¬ 
titions.  Each  class  contains  only  droplets  within  a  certain  range  of  sizes,  which  are  approximated 
as  one  average  size.  The  continuity  equation  for  droplets  of  class  k  is  very  similar  to  that  for  the 
gaseous  species  and  is  solved  in  a  fully-coupled  approach  with  the  fluid  dynamics.  The  class-fc  droplet 
continuity  equation  is  implemented  in  GASP  v4  as: 

Here  pk  is  the  density  of  droplet  class  k  in  the  two-phase  fluid.  The  droplet  source  term,  p*,  is 
modeled  as  three  separate  components 

Pk  =  Pk,  1  +  Pk, 2  +  Pk,3  j  (35) 

which  represent  nucleation,  growth,  and  transport  of  droplets  between  the  classes  respectively. 


6.1  Nucleation 

The  nucleation  source  term,  p ^ ,  is  given  as 

.  _  f  PlV*  J  :  for  class  k  where  i/2  <  r*  <  rfc+1/2 

^kl  ~  \  0  :  all  other  classes ,  ^ 

where  pl  is  the  density  of  the  liquid  water,  V*  (=  4/37T r*3)  is  the  volume  of  a  nucleating  droplet,  J 
is  the  nucleation  rate  per  unit  time  and  volume,  and  r*  is  the  critical  radius  of  a  nucleating  droplet. 
The  nucleation  rate  is  given  as 


J  = 


m 


-3/2  Pv 


h2o 


PL 


exp 


where  ctq  is  the  surface  tension  for  a  plane  surface,  computed  as 


/  0.0761  +  0.000155  (273.15  -  T) :  T<  764K 
\  0.0:  T  >  764  if , 


(37) 
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rriH2o  is  the  mass  of  a  water  molecule,  pv  is  the  density  of  the  water  vapor,  k  is  Boltzmann’s  constant 
and  AG  is  given  as 


A  G  = 


mHzO 
pi  InskT 


) 


2 


<70 


(38) 


Here,  s  (=  pv/pSo )  is  the  supersaturation  ratio  where  pv  is  the  partial  pressure  of  the  water  vapor 
and  pso  is  the  saturation  pressure  of  a  plane  liquid  surface  given  as 

ps0  =  exp  (-1187.5658  +  529.80132  In T  -78.600891  In 2T  +  3.9393576  ln3T)  .  (39) 


The  critical  radius  is  given  by 

2  op 

PlRvT  In  (s)  ’ 

where  Rv  is  the  gas  constant  of  the  water  vapor. 


(40) 


6.2  Growth 

Following  the  onset  of  nucleation,  the  radii  of  existing  droplets  increase  due  to  further  condensation 
and  evaporation.  The  source  term  for  droplet  class  k  due  to  growth  is  given  as 


Pk,  2  —  71  k 


(41) 


where  rik{~  pk/PLVk )  is  the  number  density  of  class-fc  droplets.  The  rate  of  change  in  droplet  radius 
in  Eq  (41)  is  modeled  as 

dr  e  (Pv- Psr )  /42\ 

dt  PL  (2ttRvT)1/2  ’ 

where  6  is  the  condensation  coefficient  given  as 


f  0.5  :  T  >  270  K 

6  =  1  1.0  -  0.0125  (T  -  230)  :  230  K  <  T  <  270 K  (43) 

I  1.0  :  T  <  230  if , 


and  psr  is  the  saturation  pressure  of  a  droplet  of  radius  r,  given  as 


Psr  =  Pso  exp 


( 


2a  \ 
pL  RvTr) 


(44) 


Here  cr(=  £ao)  is  the  surface  tension  of  the  droplet  and  £  is  the  coefficient  of  surface  tension.  Note 
that  in  our  calculations  £  =  1.0. 


6.3  Transport 

Droplet  growth  also  produces  a  transport  of  droplets  ( i.e a  droplet  flux)  across  the  partition 
boundaries.  For  example,  when  droplets  at  the  k  +  1/2  partition  boundary  grow  they  leave  class 
k  and  become  part  of  the  k  +  1  class.  In  order  to  model  the  number  density  distribution  within 
a  class,  we  introduce  a  distribution  function  of  the  number  density,  z/r,  such  that  the  number  of 
droplets  with  radii  between  r  and  r  4-  dr  is  given  by  dn  =  vT  dr.  Utilizing  the  distribution  function, 
the  source  term  for  class  k  due  to  transport  is  given  as 

PM  =  (PL  V  (|)  t_i/s  -  (PL  V  w  (I)  •  <«) 


23 


Notice  that  growth  at  the  k- 1/2  partition  boundary  gives  a  positive  contribution  to  the  source  term 
while  growth  at  the  k  + 1/2  boundary  gives  a  negative  contribution.  Utilizing  a  linear  reconstruction 
for  the  grouping  (pL  V  iv)fe±1/2  and  further  assuming  that  vT  is  a  piecewise  linear  function  yields 


(PlVut) 


_  1 

fc-i/2  -  Ark-i  +  A rk 


Pk- 1 
Arfc_i 


+ 


Arjfc-x 


Pk 
A  rk 


(46) 


{pL  V  Vr)k+i/2 


1 

A  rk  +  Arfc+i 


Ar‘+1£;+Ar‘ 


Pfc+i 

Arfc+iJ  ’ 


(47) 


where  A rk  =  rfc+1/2  -  rk_1/2. 


6.4  Coupling  Issues 

The  liquid-phase  species  are  treated  in  a  nearly  identical  fashion  to  the  gas-phase  species.  The 
mixture  density  for  the  two-phase  fluid  is  determined  as  the  sum  of  the  specie  densities  for  both  the 
gas  and  liquid-phase  species  as  follows 


N.  Nk 

p  =  ^_jpi  +  '52pk-  (48) 

i=l  k= 1 

Here,  Ns  is  the  total  number  of  gaseous  species  and  Nk  is  the  total  number  of  liquid  species.  The 
mixture  pressure  is  determined  as  the  sum  of  the  partial  pressures  of  only  the  gaseous  species 


P  =  f^Pi  =  f^PiRiT  =  PRT,  (49) 

i=  1  i=l 

where  the  mixture  gas  constant,  R ,  is  given  as 


N’  ^Ri.  (50) 

P 

The  energy  for  the  liquid-phase  species  is  given  as 

ek=  f  CVLdr  +  +  C  k  =  1,  2, Nk ,  (51) 

where  C  is  the  latent  heat. 

The  liquid-phase  species’  continuity  equations  do  not  include  viscous  diffusion  terms  and  the 
liquid  species  do  not  contribute  to  the  mixture  viscosity,  mixture  thermal  conductivity,  or  effective 
diffusion  coefficients.  To  account  for  the  interphase  mass  transport,  the  chemical-production  source 
term  for  water  vapor  is  modified  to  include  the  effects  of  condensation  or  evaporation,  i.e.,  pv  — 

TT^Nk  * 

Pk • 

6.5  Power-Off  RADICL  Simulation  with  Condensation 

Two  additional  RADICL  simulations  have  been  performed  using  GASP  v4  to  study  the  effects 
of  water- vapor  condensation.  These  simulations  were  both  obtained  for  power-off  conditions  and 
include  a  “baseline”  case  without  water-vapor  condensation  and  a  “condensation-on”  case  that 
includes  condensation  effects.  Both  calculations  were  performed  using  the  same  flow  conditions  and 
grid  system  as  the  previously-described  RADICL  simulation  (see  Section  5.2). 
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The  condensation  calculation  utilized  5  water,  classes  with  average  droplet  radii  of 


n 

=  7  x  10-lom 

T2 

=  2  x  10-9  m 

rz 

=  1  x  10~8  m 

u 

=  1x10  ~7m 

r5 

=  1  x  10~6  m . 

(52) 

Condensation  occurs  in  the  diverging  section  of  the  nozzle,  where  the  rapid  decrease  in  temperature 
results  in  the  nucleation  and  subsequent  growth  of  water  droplets.  Figure  16  shows  contours  of  the 
log  of  the  number  density,  i.e.,  Log10(nfc),  for  each  of  the  droplet  classes  in  the  k  =  5  vertical  plane. 
Nucleation  begins  just  downstream  of  the  throat,  and  occurs  entirely  in  the  7.0  x  10“lom  class. 
Droplet  growth  and  transport  occurs  rapidly,  and  significant  number  densities  occur  in  the  larger 
droplet  classes.  This  process  occurs  progressively  further  down  stream  for  each  class  due  to  the 
finite  time  scale  for  droplet  growth  and  transport.  The  largest  1.0  x  10_6m  class  only  has  significant 
number  densities  further  down  stream  near  the  optical  cavity. 

Figure  17  shows  the  Log10(n&)  profile  for  each  of  the  water  droplet  classes  at  x  =  0.11578m 
(i.e.,  the  optical  axis)  in  the  k  =  5  vertical  plane.  The  vertical  axis  represents  the  y  coordinate 
and  starts  at  the  symmetry  plane  (i.e.,  y  =  0.0)  and  ends  at  the  upper  nozzle  wall.  Notice  that 
practically  no  water  droplets  form  in  the  boundary  layer  near  the  upper  wall.  The  GASP  v4 
simulation  predicted  maximum  number  densities  of  n\  ~  1019*5,  n 2  ~  1017,6,  n$  ~  1014,3,  n 4  ~ 
108'9,  and  n5  ~  102  6 particles/m3.  These  number  densities  are  very  similar  to  those  predicted  by 
Masuda  [10]  for  the  S-Coil  system. 

Figure  19  shows  gain  contours  for  the  condensation-on  case  and  the  baseline  case.  The  gain  con¬ 
tours  are  very  similar  in  the  upstream  region,  however,  as  condensation  begins  to  become  significant, 
the  gain  is  reduced  for  the  condensation-on  case.  This  fact  is  evident  by  the  diminished  extent  of  the 
high-gain  (i.e.,  red)  contours  for  the  condensation-on  case.  As  the  water  vapor  condenses,  the  latent 
heat  release  into  the  flow  raises  the  temperature.  This  increase  in  temperature  causes  a  significant 
reduction  in  the  concentration  of  excited  iodine,  I* ,  and  therfore  reduces  the  gain.  To  illustrate  this 
effect,  consider  the  pumping  reaction  which  controls  the  production  of  excited  iodine: 

/  +  02(1A)^7*  +  02(3E)  (53) 

where 

kf  -  4.697  x  1010  (m3/Kg-mol/s)  (54) 

kt  =  6.323  x  1010  exp~401-4/T  (m3/Kg-mol/s) .  (55) 

The  backward  rate  increases  with  increasing  temperature,  and  leads  to  the  reduction  in  excited 
iodine.  [Note:  In  GASP  the  backward  reactions  are  computed  using  the  LeRC  curve  fits  coupled 
with  the  minimization  of  Gibb’s  free  energy.  However,  the  backward  Arrhenius  rate  in  Eq  (55) 
represents  an  equivalent  formulation.]  Figure  18  shows  the  gain  and  temperature  profiles  at  the 
x  =  0.11578m  station  for  both  the  baseline  and  condensation-on  case.  The  gain  is  about  4%  lower 
in  the  core  flow  region  for  the  condensation-on  case  compared  to  the  baseline  case.  The  temperature 
in  the  core  flow  region  for  the  condensation-on  case  is  approximately  7 K  higher  than  for  the  baseline 
case.  Because  the  gain  is  a  measure  of  the  flow’s  photon  extraction  potential,  we  would  expect  a 
reduction  in  the  output  power  for  the  condensation-on  case. 

Initial  attemps  at  running  the  RADICL  configuration  with  condensation  yielded  much  less  liquid 
water  than  than  was  experienced  by  Masuda  [10]  for  similar  conditions.  Examination  of  the  results 
indicated  that  the  solution  could  be  greatly  affected  by  the  droplet  sizes  chosen.  Initially  we  chose 
droplet  radii  of  rx  =  1  x  10_1°  m  and  r2  =  lx  10”9  m  for  the  first  two  classes.  The  remaining  three 
classes  had  the  same  radii  as  described  in  Eq  (52).  For  this  case,  the  first  droplet  class  had  partition 
boundaries  of  5.5  x  10_11m  on  the  lower  boundary  and  5.5  x  10~10m  on  the  upper  boundary.  It 
turned  out  for  this  droplet  discretization,  most  of  the  nucleating  droplets  evaporated  before  they 
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could  grow  and  transport.  To  see  this,  we  rewrite  Eq  (42)  in  terms  of  the  supersaturation  ratio,  s, 
and  the  critical  radius,  r* 

dr  _  6>Ml-*rVr—1) 

9t~  PL(2nRvTf2  ’  (56) 

The  sign  of  drjdt  depends  on  the  size  of  the  droplet  in  relation  to  the  critical  radius  as  follows: 

r  =  r*  -4  r*  jr  -1  —  0  ->  drjdt  —  0 

r  <  r*  ;  s  >  1  ->  r*  jr  -  1  >  0  -4  <9r/<9£  <  0 

r  >  r* ;  s  >  1  -4  r*/r  —  1  <  0  -4  <9r/<9£  >  0 

Note  that  when  s  >  0  the  vapor  pressure  is  greater  than  the  saturation  pressure  and  condensation 
should  occur.  For  the  RADICL  conditions,  the  critical  radius  at  which  droplets  nucleate  was  within 
the  range  of  1.69  x  10-lom  <  r*  <  6.5  x  10-lom.  For  the  initial  droplet  discretization,  this 
range  of  critical  radii  occurred  almost  entirely  between  the  average  droplet  size  of  the  first  class 
(■ i-e .,  1  x  10_10m)  and  the  right  partition  boundary  ( i.e .,  5.5  x  10-lom).  For  this  case  r\  <  r*  and 
dri/dt  <  0.  Hence  the  droplets  that  nucleate  into  the  first  class  begin  immediately  to  evaporate. 
The  rate  at  which  the  droplets  evaporate  depends  on  the  magnitude  of  the  supersaturation  ratio. 
For  the  RADICL  conditions  under  consideration,  the  supersaturation  ratio  is  large  and  on  the 
order  of  s  =  106,  causing  rapid  evaporation  of  nucleating  droplets.  For  the  modified  set  of  droplet 
radii  described  in  Eq  (52),  nucleation  occurs  entirely  between  the  left  partition  boundary  (i.e., 
1  x  10_lom)  and  the  average  droplet  size  for  the  first  class  (i.e.,  7  x  10“10m).  For  this  case  r\  >  r* 
and  dri/dt  >  0.  Hence  the  droplets  that  nucleate  into  the  first  class  begin  immediately  to  grow. 
Increasing  the  number  of  droplet  partitions  can  help  alleviate  this  discretization  problem,  however, 
future  research  must  be  conducted  to  determine  a  logical  way  to  select  droplet  sizes  and  partition 
boundaries. 


7  Wall  Catalysis  Model 

COIL  efficiency  can  also  be  adversely  affected  by  wall-catalysis  reactions.  These  reactions  promote 
the  deactivation  of  excited  species  that  collide  with  the  nozzle  wall.  Wall  catalysis  effects  on  the 
gas-phase  chemical  composition  are  modeled  in  GASP  v4  using  a  surface  mass-balance  boundary 
condition  coupled  with  specified  wall  deactivation  and  recombination  reactions.  The  wall  mass 
balance  equation  for  the  ith  specie  is  given  as 

Nr 

PiVi-  n  =  i  =  1,2, ...,Na,  (57) 

r= 1 

where  Vi  is  the  i^/i-specie’s  diffusional  velocity  and  o;[  represents  chemical  production  of  the  ith 
specie  due  to  surface  reaction  r.  Equation  (57)  states  that  the  diffusional  mass  flux  of  specie  i 
leaving  the  wall  must  equal  the  rate  of  production  of  specie  i  through  wall  catalysis  reactions.  To 
be  consistent,  the  wall  diffusion  model  utilized  in  Eq  (57)  should  be  the  same  as  that  used  in  the 
interior  flow,  i.e.,  the  effective  diffusion  model. 


7.1  Surface  Mass  Balance 


The  surface  mass  balance  equation  has  been  implemented  within  the  framework  of  a  fixed-temperature, 
no-slip,  wall  boundary  condition.  The  no-slip,  wall  boundary  condition  utilizes  a  zero-pressure- 
gradient  condition  normal  to  the  wall  to  determine  the  wall  pressure.  As  a  result  of  this  assumption, 
the  pressure  diffusion  terms  in  Eq  (16)  can  be  neglected,  and  substitution  of  the  effective  diffusional 
mass  flux  into  Eq  (57)  yields 


Pi 
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3= 1 


r= 1 


i  =  1, 2, ...,  Ns 


(58) 
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Figure  17:  Water  droplet  number  density  profile,  Log10(nfc),  profiles  at  the  RADICL  optical  axis 
(x  =  0.11578 m). 


(a)  Gain  profile.  (b)  Temperature  profile. 


Figure  18:  Gain  (1/m)  and  temperature  (K)  profiles  at  the  RADICL  optical  axis  (x  =  0.11578m) 


r 

Reaction 

Pi  (FC) 

PI  (Mod) 

1 

O2  0  E)  +  wall  - *  O2  (x  A)  +  waN 

1.0 

0.0 

2 

O2CE)  +  wa^  ^2(3E)  +  waii 

0.0 

1.6e“2 

3 

02  0 A)  +  wall  ->  02(3E)  +  wall 

1.0 

1.3e-5 

4 

/(2 Pl/2)  +  wall  -4  /(2P3/2)  +  wall 

1.0 

1.0 

5 

I2  +  wall  I2+  wall 

1.0 

1.0 

6 

2  I(2Px/2)  +  wall  ->  I2  +  wall 

0.0 

0.0 

7 

2  /(2P3/2)  +  wall  — y  I2  T  wall 

1.0 

1.0 

Table  4:  Catalytic  wall  reactions  with  fully-catalytic  (FC)  efficiencies  as  used  by  Hishida  [8],  and 
modified  efficiencies  (Mod)  as  suggested  by  Madden  [9]. 


Utilizing  Eq  (61),  the  catalytic  production  terms  on  the  right-hand  side  of  Eq  (58)  in  Kg/m2/s  are: 


5Zwo2(3S)  -  ^o2(3E)  (voii'T,)  +7?o2(1A)) 

r= 1 
Nr 

X)wo2(i A)  =  Mo2(  1A)  (flb2(1'Z)  -  vbit'A)) 

r== 1 

X^wo2(iE)  =  (_,?02(1E)  “  ^(‘E) j 

r=l 

XlW/(2f3/2)  =  Ml(2P3/2)  (Vl(2p1/a)  ~  Wl(*P3/2)) 
r=  1 

X^w/(2P,/2)  =  ^I(.2Pi/2)  (_7?/(2^i/2)  ~rf(ip3/2)) 

r= 1 

/  i  ^  \ 

=  Mj2  ( ?7/*  +  2^(2A/2)  +  2^ PS/2) ) 

r=1  '  ' 


r=l 


The  source  terms  for  the  non-catalytic  species,  «.e.?  #e,  #20,  and  C72 ,  are  zero.  The  species’ 
source  terms  above  are  subsititued  into  their  corresponding  mass  balance  equation,  z.e.,  Eq  (58), 
and  division  through  by  the  mixture  concentration,  77*,  yields  a  right-hand  side  in  terms  of  mole 
fractions  at  the  wall. 

In  the  most  general  procedure  Eq  (58)  is  discretized  using  finite  differences  and  solved  in  a  coupled 
fashion  for  the  Ns  species’  mole  fractions  at  each  point  on  the  surface.  First-order  discretization  of 
the  mole-fraction  gradients  in  Eq  (58)  yields 


Mi 


E;=i  x\w)Mi  pi 


E  M,Dim  (xf  -  xf )  -  (*S‘’  -  *<">) (62) 


Ah/2 
* 

i  =  1,2, Ns 


where  the  superscript  (w)  indicates  evaluation  at  the  wall,  the  superscript  (1)  indicates  evaluation 
at  the  first  cell  center  off  the  wall,  and  A/i/2  is  the  distance  between  the  wall  boundary  location  and 
the  first  cell  center.  Here  we  have  also  substituted  the  relation 


Ei  -  Xi-M-i 
p  E(=1  XI  Ml  ' 
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For  fixed  temperature  and  pressure,  Eq  (62)  represents  a  Ns  x  Ns  nonlinear  system  of  equations  in 
terms  of  the  wall  mole  fractions,  The  resulting  system  of  equations  is  solved  using  Newton’s 
method  with  an  LU  decomposition  and  forward/ back  ward  substitution  that  vectorizes  over  all  points 
on  the  surface. 

An  overview  of  the  catalytic  wall  boundary  condition  is  as  follows: 

1.  The  wall  temperature  is  fixed  to  Tw. 

2.  The  wall  pressure  is  computed  from  (Vp)^  =  0. 

3.  For  a  fixed  temperature  and  pressure,  the  nonlinear  system  given  by  Eq  (62)  is  solved 

for  the  wall  species’  mole  fractions,  Starting  values  for  the  mole  fractions  are 

taken  from  the  previous  iteration’s  converged  values. 

4.  The  total  number  of  moles  at  the  wall,  r)[w\  is  computed  from  the  thermally-perfect 
gas  law,  p  =  rjtRT. 

5.  The  wall  species’  concentrations  are  determined  using  77;  =  XiVt- 

6.  The  wall  species’  densities  are  determined  using  pi  =  rfo  Mi . 

7.  The  wall  velocity  components  are  set  to  zero,  i.e.,  u,v,w  =  0. 


Initial  attempts  to  solve  Eq  (62)  revealed  that  this  system  was  ill-behaved.  When  solutions  to 
the  mass-balance  system  could  be  obtained,  there  was  no  mechanism  to  ensure  that  the  sum  of  the 
mole  fractions  at  the  wall  would  be  unity,  i.e.,  YlXi^  =  1-  To  satisfy  this  constraint,  we  introduced 
a  penalty  function  to  each  mass  balance  equation  of  the  form 

minimize  fa  (x)  =  Fi  (x)  +  rh2  (x) 
such  that  h(x)  =  0 

where  Mx)  =  ~  1 '  (64) 

The  penalty  function  associates  a  penalty  with  a  violation  of  the  constraint.  Solutions  of  the  aug¬ 
mented  system  appears  to  be  well-behaved,  satisfying  both  the  initial  system  as  well  as  the  mole- 
fraction  constraint.  Note  that  r  is  a  multiplier  which,  in  the  true  sense  of  a  penalty  function  method, 
is  systematically  increased  from  zero  to  some  significant  value  during  the  iteration  process.  In  our 
case  we  have  successfully  used  v  —  1  during  the  entire  iteration  process. 

7.2  Power-Off  RADICL  Simulations  with  Wall  Catalysis 

Two  additional  power-off  RADICL  simulations  have  been  performed  using  GASP  v4.  These  cal¬ 
culations  utilized  the  surface  catalysis  reactions  for  the  nozzle  wall  boundary  conditions.  The  first 
calculation  used  the  fully  catalytic  efficiencies  used  by  Hishida  [8]  The  second  calculation  used  the 
more  realistic  catalytic  efficiencies  suggested  by  Madden  [9].  Both  calculations  were  performed  us¬ 
ing  the  same  flow  conditions  and  grid  system  as  the  previously-described  RADICL  simulation  (see 
Section  5.2). 

Figure  20  shows  the  oxygen  species’  mass-fraction  profiles  at  the  start  of  the  constant-area  section 
upstream  of  the  injectors  (i.e.,  x  =  0.0  in  Figure  13)  in  the  k  =  5  vertical  plane.  The  vertical  axis 
represents  the  y  coordinate  and  starts  at  the  symmetry  plane  (i.e.,  y  —  0.0)  and  ends  at  the  upper 
nozzle  wall.  Results  are  shown  for  the  fully-catalytic  efficiencies,  the  modified  efficiencies,  and  the 
baseline  case.  At  this  location  there  is  no  iodine  in  the  flow  and  the  chemical  composition  is  affected 
only  by  surface  reactions  1  —  3. 

Figures  20(a)  and  20(b)  show  the  mass-fraction  profiles  for  the  ground  state  of  oxygen,  02(3^). 
Figure  20(b)  shows  a  close-up  of  the  modified-efficiency  case.  Both  catalytic  cases  predict  higher 
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near- wall  concentrations  of  02(3£)  compared  to  the  baseline  case.  The  mass-fraction  gradient  of 
02(3£)  at  the  wall  for  the  catalytic  cases  indicates  a  diffusional  flux  away  from  the  wall.  This  flux 
balances  with  the  production  of  02(3£)  through  surface  deactivation  reactions  2  and  3.  As  expected, 
the  fully-catalytic  efficiencies  have  a  much  greater  impact  on  the  near- wall  concentration  of  02(3£) 
compared  to  the  modified-efficiency  case. 

Figures  20(c)  and  20(d)  show  the  mass-fraction  profiles  of  the  02(xA)  state  of  oxygen.  Fig¬ 
ure  20(d)  shows  a  close-up  of  the  modified-efficiency  case.  Both  catalytic  cases  predict  lower  near- wall 
concentrations  of  02(*  A)  compared  to  the  baseline  case.  The  diffusional  flux  of  02(XA)  toward  the 
wall  balances  with  the  production  and  destruction  of  02(1  A)  through  surface  deactivation  reactions 
1  and  3.  Although  02(xA)  is  produced  in  surface  reaction  1,  this  process  is  outweighed  by  reaction 
3  and  the  net  effect  is  to  destroy  02(XA)  at  the  wall.  Again,  the  fully-catalytic  efficiencies  have  a 
much  greater  impact  on  the  near- wall  concentration  of  02(XA)  compared  to  the  modified-efficiency 
case. 

Figures  20(e)  and  20(f)  show  similar  mass-fraction  profiles  of  the  02(XS)  state  of  oxygen.  Both 
catalytic  cases  predict  lower  near- wall  concentrations  of  02(XS)  compared  to  the  baseline  case, 
resulting  from  the  fact  that  02(1E)  is  destroyed  through  surface  deactivation  reactions  1  and  2. 

Figure  21  shows  mass  fraction  profiles  for  each  of  the  iodine  species  at  the  optical  axis  (i.e., 
x  —  0.11578m  in  Figure  13)  in  the  k  =  5  vertical  plane.  Again,  results  are  shown  for  the  fully- 
catalytic  efficiencies,  the  modified  efficiencies,  and  the  baseline  case. 

Figures  21(a)  shows  the  mass-fraction  profiles  for  the  ground  state  of  diatomic  iodine,  I2.  Both 
catalytic  cases  predict  higher  near-wall  concentrations  of  I2  compared  to  the  baseline  case.  The 
mass-fraction  gradient  of  72  at  the  wall  for  the  catalytic  cases  indicates  a  diffusional  flux  away  from 
the  wail.  This  flux  balances  with  the  production  of  72  through  surface  deactivation  reaction  5  and 
recombination  reactions  6  and  7.  The  core-flow  concentration  of  72  is  significantly  higher  for  the 
fully-catalytic  case  compared  to  the  baseline  and  modified-efficiency  calculations.  In  this  case  the 
very  strong  upstream  catalytic  effects  have  propagated  into  the  core  flow  and  altered  the  gas-phase 
chemical  mechanisms.  The  lower  upstream  concentrations  of  02(XE)  and  02(x  A)  lead  to  a  reduction 
in  the  dissociation  of  72  that  occurs  through  gas-phase  reactions  7  and  11.  The  net  effect  is  to  have 
more  J2  in  the  core  flow.  Differences  between  the  modified-efficiency  and  the  baseline  case  are 
confined  to  the  near- wall  region. 

Figure  21(b)  shows  the  mass-fraction  profiles  for  the  excited  state  of  diatomic  iodine,  7£  •  Both 
catalytic  cases  predict  lower  near-wall  concentrations  of  7£  compared  to  the  baseline  case.  The 
diffusional  flux  of  7£  toward  the  wall  balances  with  the  destruction  of  7 J  through  surface  deactivation 
reaction  5.  Again,  the  core-flow  concentration  of  1%  significantly  higher  for  the  fully-catalytic 
case  compared  to  the  baseline  and  modified-efficiency  calculations.  The  reason  for  the  increased 
concentration  of  7|  is  the  reduction  in  amount  of  J2  dissociation  (as  described  above)  and  therefore 
the  overall  increase  in  the  amount  of  72  available  to  become  excited.  Again,  differences  between  the 
modified-efficiency  and  the  baseline  case  are  confined  to  the  near-wall  region. 

Figure  21(c)  shows  the  mass-fraction  profiles  for  the  ground  state  of  monatomic  iodine,  7(2P3/2)« 
Both  catalytic  cases  predict  lower  near- wall  concentrations  of  7(2P3/2)  compared  to  the  baseline  case. 
The  mass-fraction  gradient  of  I(2P3/2)  at  the  wall  for  the  catalytic  cases  indicates  a  diffusional  flux 
toward  the  wall.  This  flux  balances  with  the  production  of  7(2P3/2)  through  surface  deactivation 
reaction  4  and  the  destruction  of  7(2P3/2)  through  surface  recombination  reaction  7.  Note  that 
surface  reaction  7  dominates  reaction  4  so  that  the  net  effect  is  to  destroy  7(2P3/2)  at  the  wall.  The 
core-flow  concentration  of  7(2P3/2)  is  significantly  lower  for  the  fully-catalytic  case  compared  to  the 
baseline  and  modified-efficiency  calculations.  The  reason  for  the  reduced  concentration  of  7(2P3/2) 
is  again  the  reduction  in  amount  of  72  dissociation.  Differences  between  the  modified-efficiency  and 
the  baseline  case  are  confined  to  the  near-wall  region. 

Figure  21(d)  shows  the  mass-fraction  profiles  for  the  excited  state  of  monatomic  iodine,  7(2P!/2). 
Both  catalytic  cases  predict  lower  near-wall  concentrations  of  7(2Px/2)  compared  to  the  baseline  case. 
The  diffusional  flux  of  I(2Pi/2)  toward  the  wall  balances  with  the  destruction  of  7(2Pl/2)  through 
surface  deactivation  reaction  4  and  surface  recombination  reaction  6.  The  core-flow  concentration  of 
7(2Pi/2)  is  significantly  lower  for  the  fully-catalytic  case  because  of  the  reduction  in  the  amount  of 
72  dissociation.  Again,  differences  between  the  modified-efficiency  and  the  baseline  case  are  confined 
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Figure  22:  Gain  distribution  (1/m)  at  the  RADICL  optical  axis  ( x  —  0.11578m). 


more  to  the  near-wall  region. 

Figure  22  shows  the  gain  profiles  at  the  optical  axis  (x  =  0.11578m)  in  the  k  =  5  vertical  plane. 
Both  catalytic  cases  predict  lower  values  of  gain  compared  to  the  baseline  case.  The  differences 
between  the  modified-efficiency  and  the  baseline  case  is  confined  to  the  near-wall  region  while  the 
fully-catalytic  case  effects  the  gain  in  the  core  flow.  Because  the  gain  is  a  measure  of  the  flow’s 
photon  extraction  potential,  we  would  expect  a  reduction  in  the  output  power  for  the  wall-catalysis 
cases.  This  effect  will  be  much  larger  for  the  fully-catalytic  case. 

7.3  Power-On  RADICL  Simulation  with  Water- Vapor  Condensation  and 
Wall  Catalysis 

A  power-on  RADICL  simulation  has  been  performed  with  both  water-vapor  condensation  and  sur¬ 
face  catalysis.  The  simulation  considered  the  same  5  droplet  classes  as  described  in  the  previous 
power-off  condensation  case  and  the  modified  efficiencies  were  utilized  for  the  surface  catalysis  re¬ 
actions.  Nozzle  and  injector  inflow  conditions  were  the  same  as  the  previous  cases.  The  unit-cell 
approximation  was  again  employed  with  the  same  grid  topology. 

The  aperture  (and  associated  mirror  system)  was  located  between  x  =  0.08578  m  and  x  = 
0.14578  m.  The  length  of  the  aperture  was  La  =  0.06  m.  The  distance  between  the  spherical  mirror 
and  the  flat  out-coupling  mirror  is  Dm  =  0.8  m  and  the  radius  of  curvature  of  the  spherical  mirror  is 
r2  =  10.0  m.  The  coefficients  of  reflectivity  are  Ri  =  0.927  and  R2  =  0.995  for  the  flat  out-coupler 
and  spherical  mirrors  respectively.  The  gain  length  was  taken  as  Lg  =  0.254  m  (the  lateral  width 
of  the  nozzle).  Distributed  losses  in  the  gain  medium  and  diffractive  losses  at  the  aperture  were 
neglected.  The  polar  mesh  size  for  the  ray-tracing  algorithm  was  101  x  41  and  the  laser  optics 
resonator  model  was  called  on  every  third  flow-solver  iteration.  A  baseline  power-on  solution  has 
been  run  for  this  case  without  condensation  and  surface  catalysis  and  results  have  been  presented 
in  Eppard  [4]  (see  Section  5.2). 

The  laser  power  output  predicted  for  the  baseline  simulation  was  P  =  7.524  KW.  With  conden¬ 
sation  and  surface  catalysis  the  simulated  power  output  was  P  =  6.538  KW\  a  13%  reduction  in 
power  due  to  the  combined  effects  of  condensation  and  surface  catalysis.  The  experimentally  mea¬ 
sured  value  of  power  was  P  =  5.8  KW.  Our  results  are  in  reasonable  agreement  with  experiment, 
since  no  optical  loss  terms  were  included.  Similar  studies  ( e.g Buggeln  [1])  have  used  non-zero 
values  for  mirror  scattering  and  diffractive  losses  to  give  corrected  output  powers  that  are  up  to  20% 
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lower  than  the  zero-loss  power.  Assuming  similar  behavior  for  this  case,  our  power  results  are  in 
good  agreement  with  the  measured  value.  The  results  of  this  study  were  presented  in  Eppard  [3]. 


8  COIL  Sensitivity  Analysis 

The  sensitivity  implementation  in  SENSE  is  fundamentally  based  on  implicit-differentiation  applied 
to  the  boundary- value  problem  describing  the  fluid  mechanics.  That  is,  one  envisions  an  implicit 
equation 

n  Q,/?)  =  0,  (65) 

where  Q  denotes  the  distributed  dependent  variables  ( e.g .  density,  momentum,  energy)  and  p 
denotes  a  design  variable.  In  common  terminology  Q  is  known  as  the  state  and  /?  as  the  design 
parameter.  For  fixed  /?,  Eq  (65)  is  solved  for  Q.  This  defines  a  map 


which  associates  a  flow  solution  with  the  specified  design  parameter (s).  The  sensitivity  is  the  deriva¬ 
tive  of  this  map;  it  provides  a  linear  approximation  for  how  the  flow  solution  will  change  under  a 
small  change  in  the  design  parameter.  One  proceeds  by  formally  differentiating  Eq  (65)  to  produce 


dlld  Q  dll 
5Q  dp  +  dp 


(66) 


The  SENSE  code  implements  a  numerical  approximation  to  Eq  (66)  for  situations  wherein  H  models 
Reynolds-Averaged  Navier-Stokes  flows  with  finite-rate  chemistry. 

The  integral  equations  for  the  fluid  dynamic  system  are  written  as 


it  III qdv + £ (F(Q)  • dA = £  (f"(q)  • fi) dA + HI s dy>  ( 6 7) 


where  Q  =  Q (x,y,z,t;P)  represents  the  vector  of  state  variables  resulting  from  conservation  of 
mass,  momentum  and  energy.  The  surface  integrals  represent  the  inviscid  and  viscous  fluxes  (F  and 
F„).  Changes  in  chemical  composition  are  governed  by  the  species  production  terms  in  the  source 
term  S.  Formal  differentiation  of  Eq  (67)  with  respect  to  the  generic  design  parameter,  /?,  results  in 
the  sensitivity  equations.  These  equations  are  linear  in  the  sensitivities.  Presented  in  integral  form, 
they  are 

i  III Ql  dV + f A  (F,(Q)  -h)dA  =  fA  (F-(Q)  'A)dA+ IIIs' (68) 

where  Q'  is  the  flow  sensitivity,  i.e.,  dQ/d/3.  The  inviscid  and  viscous  flux  sensitivities,  F'  and  F'v, 
and  the  source  term  sensitivity,  S',  are  formally  given  as 

9F  9Q  _  3F.  dQ  #  dS  dQ 

SQ  d/3  ’  v  9Q  d/?  ’  3Q  d/T  V  ; 


SENSE  solves  Eq  (68)  using  an  upwind  characteristic-based  formulation. 


8.1  Sensitivities  for  the  Effective  Diffusion  Model 

As  a  first  step  toward  tailoring  SENSE  for  COIL-laser  applications,  AeroSoft  has  extended  the 
viscous  sensitivity  formulation  to  include  the  effective  diffusion  model.  This  task  required  modifying 
the  species  diffusive  flux  sensitivity  terms  which  are  given  as 

F'v  =  --^{psVs-n)  ,  (70) 
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where  the  diffusion  velocity  is  determined  by  Eqs  (16)  (19)  and  (20).  The  sensitivity  formulation 
proceeds  by  formally  differentiating  Eqs  (19)  and  (20)  with  respect  to  the  design  variable,  yielding 


(p,V?)W  Gf-^gGJ 


,  (71) 


(**)'  =  -7.'  ^1  a  lor  -  i  (<4-  */)  Eg?  -  ^'XXj  f 


p/  Vp 


g?-tEg; 


VP„'4.  VP'^ 


where 


The  individual  sensitivity  terms  for  the  species  concentration,  the  mixture  concentration,  and  the 
species  mole  fraction  are  given  as 


V  =  fv  v'=ifA_. 

7t  2_,7S>  7t  \MS  ‘ 


The  sensitivities  of  G*,  G§,  and  £>Sfn  are 


G*'  =  Ms  D'sm  VXs  +  MsDsmVx's 


g (x.-^)  +*»c. 


D'  = 
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Finally,  closure  of  the  effective  diffusion  sensitivity  terms  is  given  by  sensitivities  of  the  binary 
diffusion  coefficient 

w  _*vV3  T'  %  P'\  rn\ 


where 


v'ij  Vij  V  2  t  a 


0.06606  0.57901  1  *, 

ji*1.08742  "I"  j.*1-84910  J  > 


m„  V 


The  sensitivity  relations  given  by  Eqs  (71)  -  (80)  along  with  the  corresponding  flux-sensitivity 
Jacobian  terms  have  been  added  to  a  COIL  version  of  the  SENSE  software.  A  simplified  helium- 
iodine  chemistry  model  has  also  been  added  to  the  SENSE  package. 
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8.2  Sensitivities  of  Helium-Iodine  Injection 

The  sensitivity  analysis  has  been  applied  to  the  helium-iodine  injection  case  described  in  Section  3.1. 
For  this  case  the  design  variable  was  chosen  to  be  the  jet  iodine  density,  i.e .,  /?  =  pi2Jet .  Fig¬ 
ures  23(a)  and  23(b)  show  the  iodine  mass  fractions  from  the  GASP  calculation  and  the  iodine 
mass-fraction  sensitivities,  (pi2/pY  from  SENSE.  The  iodine  mass-fraction  sensitivities  at  the  in¬ 
jector  face  are  imposed  by  the  jet  inflow  densities  and  density  sensitivities  through  the  following 
relation 


The  jet  helium  and  mixture  densities  are  p/2  =  0.034737  Kg/m ?  and  p  =  .082217  Kg/mz .  The  jet 
helium  and  mixture  density  sensitivities  are  unity  for  this  case.  Under  these  conditions  Eq  (81)  gives 
a  value  of  (pj2/p)'  =  7.024;  the  free-stream  value  observed  in  Fig  23(b).  The  iodine  mass-fraction 
sensitivities  are  observed  to  have  values  larger  than  {pi2/p)'  —  7.024  near  the  injector  side-wall 
regions.  This  indicates  that  the  pressure-diffusion  effects  become  more  pronounced  as  the  iodine 
density  in  the  jet  increases.  To  see  that  this  is  plausible,  we  analyze  the  leading  pressure-gradient 
term  in  Eq  (15),  (xi2  -  pj2/p).  For  a  5%  increase  in  the  design  variable  the  absolute  value  of  this 
term  increases  from  \(xi2  ~  Ph/p) I  —  0.387750  to  a  value  of  \{xi2  ~~  Ph/pY  ~  0.397994  (a  2.64% 
increase).  As  a  result  the  pressure  diffusion  effects  should  become  more  pronounced.  Figure  24(a) 
shows  iodine  mass  fraction  sensitivities  along  the  i  =  10  grid  line  for  both  SENSE  and  a  central 
difference  approximation  using  GASP  and  a  5%  change  in  design  variable.  The  two  profiles  agree 
very  well  and  both  predict  an  increase  in  the  iodine  mass  fraction  near  injector  side  wall.  Figure  24(b) 
shows  a  close  up  of  iodine  mass  fraction  profiles  for  the  baseline  solution,  a  GASP  solution  for  a 
5%  increase  in  /?,  and  a  Taylor’s  series  projection  using  the  baseline  flow  and  the  sensitivities  from 
SENSE.  The  projected  solution  agrees  very  well  with  the  “exact”,  GASP  solution,  except  for  the 
slight  lag  near  the  start  of  the  large  pressure-gradient  region. 

8.3  Sensitivities  of  Transverse  Helium-Iodine  Injection 

The  sensitivity  analysis  has  also  been  applied  to  the  transverse  helium-iodine  injection  case  described 
in  Section  3.2.  For  this  case  the  design  variable  was  chosen  to  be  the  jet  velocity,  i.e.,  fi  =  V^et . 
Figure  25(a)  shows  the  iodine  mass  fractions  from  the  GASP  calculation  and  Fig  25(b)  shows  the 
iodine  mass-fraction  and  velocity-vector  sensitivities.  The  velocity- vector  sensitivities  indicate  an 
increase  in  the  jet  velocity  as  well  as  deeper  penetration  of  the  helium-iodine  mixture  into  the 
primary  flow.  The  iodine  mass-fraction  sensitivities  also  indicate  that  the  jet  flow  has  penetrated 
further  into  the  primary  flow.  The  region  of  positive  (and  negative)  sensitivities  just  above  (and 
below)  the  jet  indicate  that  the  jet  has  shifted  upward.  The  portion  of  the  recirculation  region  with 
low  iodine  mass  fractions  has  also  shifted  upward,  as  indicated  by  the  positive  sensitivities  in  this 
area. 


8.4  Simple  Power-On  Test  Case 

Power-on  flow  simulations  have  been  performed  for  a  simplified  geometry,  but  with  flow  conditions 
and  optical  parameters  consistent  with  the  RADICL  laser.  This  test  problem  was  used  to  debug 
and  demonstrate  the  flow  coupling  with  the  laser  optics  resonator  model.  The  same  problem  will  be 
used  to  test  the  COIL  sensitivity  package  and  is  described  below. 

Figure  26(a)  shows  the  grid  topology  for  the  151  x  51  flow-solver  mesh  and  the  101  x  51  polar 
optics  mesh.  The  number  of  grid  points  shown  in  this  figure  has  been  reduced  to  facilitate  viewing. 
Notice  that  the  x-direction  stretching  at  the  leading  and  trailing  edges  of  the  optics  mesh  is  nearly 
identical  to  the  flow  solver  mesh.  The  start  of  the  aperture  is  at  x  =  .06  m  and  the  aperture 
length  is  La  =  .05  m.  The  aperture  height  is  Ha  =  .028  m.  The  distance  between  the  curved 
mirror  and  the  flat  out-coupling  mirror  is  Dm  =  0.8  m  and  the  radius  of  curvature  of  the  curved 
mirror  is  =  10.0  m.  The  coefficients  of  reflectivity  for  the  mirrors  are  R\  —  0.90  and  R 2  =  0.99 
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(a)  Iodine  mass  fractions. 


(b)  Iodine  mass  fraction  sensitivities. 


Figure  23:  Sensitivities  for  Helium-iodine  injection. 


for  the  out-coupler  and  curved  mirrors  respectively.  The  gain  length  was  taken  as  Lg  =  0.25  m. 
Distributed  losses  in  the  gain  medium  and  diffractive  losses  at  the  aperture  were  neglected.  The 
inflow  boundary  was  split  into  two  portions;  with  the  lower  stream  having  a  high  concentration  of 
I*  (high-gain  region),  and  an  upper  stream  with  very  little  I*  and  O2C  A)  (low-gain  region).  The 
low-gain  region  was  designed  to  account  for  the  region  of  low  gain  observed  near  the  upper  wall  in 
the  power-off  RADICL  simulations.  The  inflow  conditions  for  the  test  case  were  chosen  from  actual 
conditions  (upstream  of  the  lasing  cavity)  for  the  power-off  RADICL  simulation,  and  are  listed  in 
Table  3. 

Figures  26(b)  and  26(c)  show  the  gain  and  intensity  contours  for  this  simulation.  Values  for  the 
inflow  gain  in  the  high  and  low-gain  regions  are  approximately  a  =  1.065  1/m  and  a  =  0.00751  1/m 
respectively.  The  maximum  gain  occurs  in  the  high-gain  region  just  before  the  start  of  the  aperture 
and  has  a  value  of  amax  =  1.308  1/m.  Past  the  start  of  the  aperture,  the  gain  is  dramatically 
reduced  through  interaction  with  the  lasing  source  terms.  The  maximum  intensity  occurs  at  the 
aperture  leading  and  trailing  edges  with  a  value  of  Imax  =  3.6  x  109  VF/m2.  The  power  output  for 
this  case  converged  after  approximately  1000  iterations  and  was  P  =  13.1  KW. 

8.4.1  Power-OfF  Sensitivity  Analysis 

In  this  section  we  test  the  newly-implemented  COIL  chemistry  and  thermodynamics  sensitivity 
models.  We  consider  the  simple  geometry  described  above,  but  with  the  laser  turned  off.  In  this 
case,  the  design  variable  is  the  free-stream  density  of  O-^A)  in  the  lower  stream.  Figure  27(a) 
shows  the  gain  contours  along  the  lasing  cavity.  Figure  27(b)  shows  the  gain  distribution  along  the 
j  =  10  gridline  (midway  in  the  lower  stream).  The  gain  initially  increases  due  to  chemical  reactions 
and  reaches  a  peak  value  of  a  =  1.32  1/m.  Because  of  the  heat  release  of  the  chemical  reactions, 
and  because  the  flow  is  does  not  expand,  the  temperature  increases  causing  the  gain  to  decrease  as 
the  flow  continues  downstream.  In  our  sensitivity  analysis  we  seek  to  find  out  how  the  gain  changes 
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(a)  Iodine  mass  fraction  sensitivities  (station  1=10). 
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(b)  Projected  iodine  mass  fractions  for  5%  increase  in  (Close 
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Figure  24:  Sensitivities  and  Taylor  projection  for  helium-iodine  injection. 
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(a)  Iodine  mass  fractions. 
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(b)  Iodine  mass-fraction  and  velocity-vector  sensitivities. 


Figure  25:  Sensitivities  for  transverse  helium  iodine  injection. 
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(a)  Coarsened  flow-solver  mesh  and  polar  optics  mesh. 


(b)  Gain  contours  (1/m). 


(c)  Intensity  contours  (Watts/rn2). 


Figure  26:  Simple  test  case  grid  system  and  gain  (1/m)  and  intensity  (Watts /m2)  contours. 
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No. 

Reaction 

Kf(Kg  -  mole^m6,s) 

10 

11 

20 

h + 1*  -»•  4* + 1 

/2*  +  02(1A)  ->2J  +  02(3E) 
I*  +  H20  -J-  I  +  H20 

2.288  x  101U 

1.807  x  1011 

1.204  x  109 

Table  5:  Partial  reaction  list,  10-species,  20-reaction  COIL  chemistry  model 


with  an  increase  in  the  inflow  density  of  O2OA)  in  the  lower  stream.  Figure  27(c)  shows  the  gain 
sensitivity  distribution  along  the  j  =  10  gridline  for  the  SENSE  code  and  for  a  central  difference 
calculation  using  GASP  v4.  The  CSEM  results  are  in  very  good  agreement  with  central  difference. 

8.4.2  Power-OfF  Reaction  Rate  Sensitivity  Analysis 

In  CFD  simulations  involving  finite-rate  chemistry,  there  always  exists  some  question  concerning  the 
accuracy  of  the  reaction-rate  coefficients.  Small  variations  in  these  coefficients  can  translate  into 
significant  changes  in  the  thermo-chemical  state  of  the  flow,  affecting  distributions  of  temperature, 
pressure  and  other  flow  variables.  In  the  COIL  laser,  chemical  reactions  affect  the  lasing-cavity 
temperature,  which  in  turn  affects  the  small-signal  gain  and  lasing  efficiency.  Knowing  a  priori  which 
reaction-rate  coefficients  have  the  largest  effect  on  laser  performance  will  allow  chemists  to  prioritize 
and  possibly  narrow  their  refinement  efforts.  To  study  this  effect,  AeroSoft  has  implemented  into 
SENSE  the  capability  to  compute  flow-field  sensitivities  to  the  Arrhenius  reaction-rate  coefficients. 
The  Arrhenius  rates  are  given  as 

Kf  =  cT71e~6/kT .  (82) 

In  the  new  formulation  any  of  the  Arrhenius  coefficients  can  be  selected  as  the  design  variable.  For 
instance,  if  the  design  variable  were  chosen  to  be  the  Arrhenius  coefficient  c,  the  source  term,  S', 
would  include  a  term  like 


=  (r]+  — )  ^Lt'  +  -Kf.  (83) 

d/3  V  kT)  T  c  f 

In  this  case  the  second  term  on  the  right-hand  side  of  Eq  (83)  would  drive  the  sensitivity  problem 
and  solutions  to  this  system  of  equations  will  yield  flow-field  sensitivities  to  the  Arrhenius  coefficient 
c. 

An  analysis  of  each  of  the  reactions  in  the  COIL  model  for  the  simple  test  geometry  shows  that 
the  c  coefficient  for  reactions  10,  11,  and  20  have  the  largest  affect  on  the  gain.  The  corresponding 
reactions  are  listed  in  Table  5.  Figure  28  shows  the  gain  sensitivity  along  the  j  =  10  gridline  for 
each  of  these  reactions.  Again  the  CSEM  and  central  difference  sensitivities  are  in  good  agreement. 

8.4.3  Power-On  Sensitivity  Analysis  (Constant  Intensity) 

The  CFD /LOR  coupling  occurs  through  laser  power  extraction  terms.  Power  extraction  in  the  lasing 
cavity  manifests  itself  in  the  form  of  source  terms  in  the  iodine  species  continuity  equations,  and 
in  the  global  energy  equation.  The  species  continuity  equations  are  affected  because  components  at 
different  energy  states  are  represented  as  separate  species.  In  a  COIL  laser,  iodine  in  the  excited, 
I(2P 1/2)5  energy  state  is  stimulated  to  emit  photons  (i.e.,  stimulated  emission).  Upon  emission  of 
the  photon,  the  iodine  atom  assumes  the  lower,  I  (^3/2) ,  energy  state.  As  a  result,  power  extraction 
affects  the  concentration  of  these  two  states.  Power  extraction  also  accounts  for  a  net  energy  loss 
(through  emission  of  photons)  to  the  flow.  The  modified  species  continuity  equation  is  given  by 
Eq  (23),  where  the  last  term  on  the  right-hand  side  of  the  power-extraction  source  term.  The  source 
term  is  composed  of  the  gain,  a,  the  two-way  average  intensity,  /,  the  molecular  weight,  M5,  and  the 
energy  of  a  photon,  hv  =  90956  J/g  -  mole .  The  plus  sign  for  the  lasing  source  term  in  Eq  (23)  is 
used  for  the  ground  state  of  iodine  and  the  minus  sign  is  used  for  the  excited  state.  No  other  species 
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Figure  28:  Gain  sensitivity,  j  =  10  gridline 


in  the  model  has  a  source  term  contribution  from  power  extraction.  The  global  energy  equation  has 
a  similar  source  term,  — a  J,  to  account  for  the  energy  removed  from  the  flow. 

Similarly,  the  power  extraction  sensitivity  source  terms  take  the  form 

/ aiMs\ '  _  a' I Ms  ^  aI'Ms 
y  hv  J  hv  hv 

In  this  section  we  test  a  portion  of  the  new  COIL  source  term  sensitivities  for  the  simple  geometry 
described  above.  Here  we  assume  the  two-way  intensity  in  the  lasing  cavity  has  a  constant  value  of 
I  =  1.5  x  108  Watts/m2.  In  this  case,  the  second  term  on  the  right-hand-side  of  Eq  (84)  is  identically 
zero. 

Again,  the  design  variable  is  the  free-stream  density  of  O2O  A)  in  the  lower  stream.  Figure  29(a) 
shows  the  gain  contours  along  the  lasing  cavity.  Figure  29(b)  shows  the  gain  distribution  along  the 
j  —  10  gridline.  Past  the  start  of  the  aperture,  the  gain  is  dramatically  reduced  through  interaction 
with  the  lasing  source  terms.  Figure  29(c)  shows  the  gain  sensitivity  distribution  along  the  j  =  10 
gridline  for  the  SENSE  code  with  COIL  source  terms,  and  for  a  central  difference  calculation  using 
GASP  v4.  The  CSEM  results  are  in  very  good  agreement  with  central  difference. 


9  Coupled  Sensitivity  Analysis 

We  now  wish  to  extend  the  sensitivity  analysis  of  the  previous  sections  to  fully  include  the  coupled 
CFD/optics  model.  Specifically,  we  seek  a  solution  methodology  for  the  linear  problem  of  Eq  (66) 
when  the  residual,  i.e.,  Eq  (65),  is  a  coupled  system  involving  several  disciplines. 

We  begin  with  the  implicitly-coupled  system  and  generic  design  variable  /3  as  follows 

nf(uf,up,l3)  =  0  (85) 

Tlp(uf,up,P)  =  0.  (86) 

Here  u /  variable  describes  the  fluid  state  and  the  uv  variable  describes  laser-power  state. 
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(a)  Gain  contours  (1  /m). 
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(b)  Gain  profile  (1/m),  j  =  10  gridline. 
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(c)  Gain  sensitivity  profiles,  j  =  10  gridline. 


Figure  29:  Simple  test  case  gain  and  gain  sensitivity  for  a  constant  intensity  of  J  =  1.5  x 
10  8Watts/m2. 
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The  corresponding  coupled  sensitivity  equation  is  described  by  the  system 
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/  d  Uf  > 
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\  op  / 

One  approach  to  the  coupled  sensitivity  problem  builds  on  the  existing  SENSE  software  by  employing 
a  Jacobi  or  Gauss-Seidel  iterative  procedure  to  the  coupled  system  (87). 

An  alternative  formulation  begins  with  an  explicit  representation  of  the  laser  field  replacing  (86) 
with 


«,  =  «(*/,/?)  (88) 

The  laser  optics  resonator  model  can  be  viewed  as  a  method  for  (approximately)  evaluating  the 
function  %.  If  this  explicit  representation  is  substituted  into  (136)  we  obtain 

Kf(uf,H(uf,p),p)  =  0.  (89) 


This  result  can,  in  turn,  be  implicitly  differentiated  to  yield  the  coupled  flow-sensitivity  equation 

'<9  7 Zf  +  d  TZf  dU 
_duf  dup  duf 

Thus  far,  we  have  provided  two  forms  for  the  sensitivity  calculation:  Eq  (87)  based  on  a  general 
implicitly-coupled  formulation;  and,  Eq  (90)  for  a  partially-explicit  formulation.  The  structure 
suggested  by  the  ray- tracing  algorithm  is  an  admixture  of  these. 

Conceptually  there  are  several  distinct  intensity  fields  represented  in  the  optics  model. 

1.  an  initial  intensity  field  represented  by  J0, 

2.  a  return  intensity  field  represented  by  Ir ,  and 

3.  a  two-way  intensity  field  represented  I. 

The  purpose  of  the  optics  model  is  the  evaluation  of  the  two-way  intensity  field,  which  represents  the 
coupling  to  the  flow  code.  The  initial  and  return  intensities  are,  in  a  sense,  internal  field  variables 
used  in  the  calculation.  At  convergence,  these  two  internal  fields  have  the  same  value.  Specifically, 
the  algorithm  computes  a  return  intensity  field  and  ultimately  a  two-way  intensity  field  based  on 
the  current  gain  field  and  an  initial  intensity  field.  The  latter  is  computed  by  a  relaxation  method 
from  the  previous  invocation  of  the  routine.  Thus,  we  suggest  the  abstract  description 

u;+1=n(uf,u;,p).  (9i) 

At  convergence  u”+1  =  As  a  result,  we  use  the  following  semi-explicit  form  for  Eq  (86) 
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(90) 


Tlp(uf,up^ 3)  =up-  H(uf,Up,P). 

With  this  structure  the  coupled  sensitivity  Eq  (87)  can  be  written  as: 
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(92) 


(93) 


where  Xp  is  the  identity  operator.  Equation  (93)  is  the  approach  that  we  have  taken  in  the  COIL 
sensitivity  studies. 

The  present  study  assumes  that  the  design  parameter  of  interest  (generically,  /?)  does  not  explic¬ 
itly  appear  in  the  power-extraction  module.  As  a  practical  matter  this  rules  out  cases  where,  for 
example,  the  geometry  in  the  optical  cavity  depends  on  /?.  Mathematically,  this  assumption  means 
that  the  function  H  in  Eq  (92)  does  not  depend  explicitly  on  j3.  Our  plan  is  to  use  an  iterative 
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approach  to  solve  the  coupled  linear  system  (93).  For  example,  a  Gauss-Seidel  type  implementation 
would  be  of  the  form 
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9  Ur) 


nn 


-  cn+ 1  |  ^  ^  Q7i 

+  du: , 
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(95) 


where  we  have  used  5/  and  Sp  to  denote  the  flow  and  laser-intensity  state  sensitivities,  respectively. 
That  is, 


Sp 


d  Uf 

~d J' 

d  up 
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and 


Equation  (94)  is  solved  using  a  modification  of  AeroSoft’s  SENSE  code  that  includes  the  term 
This  accounts  for  the  effect  of  variations  of  the  laser  intensity  on  the  source  terms  in  the 
flow  equations.  The  term  is  calculated  pointwise  throughout  the  flow  field  using  the  most  recently 
calculated  intensity  sensitivity,  which  is  communicated  to  the  flow-sensitivity  code  in  the  same  way 
the  intensity  field  is  communicated  to  the  flow-solver. 

The  intensity  sensitivity,  Sp ,  is  updated  as  shown  in  Eq  (95).  The  procedure  for  implementing 
this  calculation  is  realized  in  the  new  sensitivity  module  based  in  the  geometric  ray-tracing  algorithm. 
At  an  abstract  level  the  dependence  of  the  map  H  on  the  flow  field  is  naturally  decomposed  into 
two  stages:  a  map  from  the  flow  field  to  the  gain  field;  followed  by  a  map  from  the  gain  field  to  the 
intensity  field.  In  symbols 


H(uf,up)  =  M(Af{uf),up), 


(96) 


where  V  maps  the  flow-field,  pointwise,  to  the  associated  gain  field,  and  M  maps  the  gain  field, 
intensity  field  pair  to  the  associated  intensity  field.  The  required  Frechet  derivative  is  calculated  via 
the  chain-rule 


on ^dMd u 

duf  dug  duf  ’ 


(97) 


where  ug  represents  the  gain  field. 

The  N  map,  operating  pointwise,  is  an  ordinary  function  of  the  flow-field  data.  Flow-field  sensi¬ 
tivity  information,  such  as  the  sensitivity  of  the  various  number  densities  to  the  design  parameter,  is 
available  from  the  the  modified  SENSE  code.  Given  the  functional  form  for  the  gain  and  the  chain 
rule  we  can  explicitly  calculate  the  gain-field  sensitivity.  The  intermediate  matrix  in  this  chain-rule  is 
our  calculation  of  the  term  As  in  the  coupled  flow-power  calculation  various  grid  systems  must 
be  reconciled.  Eventually,  we  compute  an  array  of  gain  sensitivities  that  represents  the  sensitivity 
of  the  gain  field  to  changes  in  the  design  parameter. 

The  function  M  maps  the  gain  field  and  initial  intensity  field  to  a  new  intensity  field  as  imple¬ 
mented  in  the  ray-tracing  code.  This  is  represented  symbolically  in  Eq  (92).  Our  initial  implementa¬ 
tion  of  the  Gauss-Seidel  iteration  for  the  sensitivity  calculation  will  follow  the  approach  used  in  the 
coupled  flow  solution.  Specifically,  the  up  variable  in  Eq  (91)  is  identified  with  the  initial  intensity 
field,  which  at  convergence,  agrees  with  the  return  intensity  field.  The  symbolic  notation  given  by 
Eq  (91)  is  implemented  in  the  ray-tracing  algorithm. 

The  implementation  of  the  second  Gauss-Seidel  update,  i.e .,  Eq  (95),  proceeds  as  follows.  The 
sensitivity  of  the  initial/  return  intensity  field  will  be  computed  as 

Sp)  +  I0(J)(R1R2)m  jz  ( es )  §  (1  -  <5P)  •  (98) 

p  p  p=i  p  p=i 
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The  /3-parameter  enters  here  through  the  gain  field  variable  ap  which  is  computed  from  the  stored 
gain  field.  The  required  expression  is 


0  (e*) 
d/3 


=  2Lges'£\ 

P=  1 


d  ap 
dP  ‘ 


(99) 


Since  ap  is  an  average  of  the  gain  along  a  particular  ray  through  the  gain  field,  it  is  a  linear  function  of 
the  gain  field  array.  Its  /3-derivative  is  this  same  linear  function  of  the  entries  in  the  gain-sensitivity 
array.  The  ^-derivative  of  the  two-way  intensity  field  is  then  computed.  The  two-way  intensity 
sensitivity  is  communicated  back  to  the  flow  sensitivity  code  SENSE  in  the  same  way  that  the 
two-way  intensity  field  is  communicated  to  the  flow  solver. 


9.1  Power-On  Sensitivity  Analysis 

In  order  to  test  this  procedure,  we  will  run  the  simple  test  problem  with  the  newly-developed  COIL 
sensitivity  package  that  includes  the  optics  resonator  sensitivities.  In  this  case  the  laser  is  turned 
on  and  the  intensity  sensitivities  are  utilized  in  the  source  terms.  Again,  the  design  variable  is 
the  free-stream  density  of  02(1A)  in  the  lower  stream.  Figure  30(a)  shows  the  gain  distribution 
along  the  j  =  10  gridline.  Past  the  start  of  the  aperture,  the  gain  is  dramatically  reduced  through 
interaction  with  the  lasing  source  terms.  Figure  30(b)  shows  the  gain  sensitivity  distribution  along 
the  j  =  10  gridline  for  the  SENSE  code  with  the  coupled  optics  sensitivity  module,  and  for  a  central 
difference  calculation  using  GASP  v4.  The  CSEM  results  are  in  very  good  agreement  with  central 
difference.  The  power  sensitivity  predicted  by  the  CSEM  method  was  Pl  =  1.01948  x  107  while 
that  predicted  by  the  central  difference  procedure  was  P '  =  1.01285  x  107.  This  represents  an  error 
of  approximately  0.65%.  Figure  31  shows  the  power  sensitivity  convergence  history  for  the  global 
scheme. 


10  Paraxial  Wave  Model 

Commonly,  laser  optics  are  modeled  by  a  ray-tracing  process,  wherein  the  path  of  a  photon  through 
the  optical  cavity  is  described  as  a  sequence  of  straight  line  segments  with  specular  reflection  at  the 
bounding  mirrors.  An  excellent  presentation  of  the  ray  tracing  methodology  is  given  in  [20],  while 
application  to  a  COIL  configuration  is  described  in  [2], 

While  it  is  convenient,  a  particle  photon  description  does  not  completely  suit  our  needs.  Certain 
physical  effects,  such  as  diffusion  (spreading  of  the  light  beam)  and  diffraction  (the  effects  of  non- 
uniform  light  speed)  are  not  well  described  by  these  models.  Moreover,  for  design  purposes  it  is 
helpful  to  know  the  sensivitity  of  the  laser  performance  to  certain  design  parameters.  Our  goals  in 
this  effort  were  to  develop  continuum  models  describing  laser  intensity  in  the  COIL  optical  cavity 
and  to  develop  tools  for  calculating  sensitivity  to  design  parameters. 

In  the  present  effort  the  optical  behavior  is  modeled  in  terms  of  a  electric  field  (wave)  propagating 
through  the  optical  cavity  (see  Figure  32).  We  model  the  electic  field  as  a  transverse  electromagnetic 
wave  (TEM),  essentially  a  perturbation  of  a  plane  wave  with  normal  (almost)  aligned  with  the  optical 
( z )  axis  of  the  optical  cavity. 

10.1  Paraxial  Wave  Optics  and  Gaussian  Beams 

As  we  are  interested  in  the  behavior  of  laser  intensity  throughout  the  optical  cavity,  we  consider  the 
physical  system  shown  in  Figure  33. 

The  schematic  shown  in  Figure  33  depicts  two  mirrors  separated  by  a  distance  D .  The  lower 
mirror  (Mi)  is  flat,  while  the  upper  mirror  (M2)  is  spherical.  Figure  33  also  illustrates  the  region 
in  which  atomic  gain  is  active.  In  the  work  presented  here,  we  specify  a  real- valued  atomic  gain 
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Figure  32:  Propagating  Wave  Approach 
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function  a(x,y,z)  as 


a{x,y,z) 


j/(z,y)  if  2  G  [al,a2]; 
JO,  otherwise. 


(100) 


In  Section  10.5  we  exploit  the  fact  that  the  gain  is  a  piecewise  constant  function  of  z. 

The  behavior  of  transverse  electromagnetic  waves  is  modeled  by  the  paraxial  wave  equation 


du  __  i 
dz  2k 


-  ik(n  —  l)]u. 


(101) 


In  equation  (101),  u(x,y,z)  denotes  the  complex  amplitude  of  the  electric  field;  the  net  atomic  gain 
is  described  by  a(x,y,z)\  and,  the  index  of  refraction  in  the  optical  cavity  is  given  by  n(x,y,z). 
Finally,  the  constant  k  is  the  wave  number:  k  >  0  denotes  a  wave  propagating  toward  increasing 
values  of  z,  and  conversely  a  negative  value  of  k  represents  a  wave  moving  toward  decreasing  values 
of  z.  The  intensity  I(x,y,z)  is  related  to  the  modulus  of  the  electric  field  according  to  I  =  \u\2. 

The  partial  differential  equation  given  by  (101)  provides  a  convenient  means  of  determining  the 
electric  field  inside  the  optical  cavity.  At  the  bounding  mirrors,  however,  a  first  principles  treatment 
requires  describing  the  evolution  of  the  electric  field  on  the  mirror  surface,  as  well  as  boundary 
conditions  relating  the  interior  and  surface  electric  fields. 

Fortunately,  in  most  quantum  electronic  applications  the  optical  beam  can  be  approximated  as 
a  transverse  electromagnetic  wave  and  represented  as  a  combination  of  Hermit  e- Gaussian  beams. 
Reflection  of  these  Hermite-Gaussian  beams  by  spherical  mirrors  (as  well  as  other  common  optical 
elements)  is  well  understood.  Excellent  presentations  of  Hermite-Gaussian  beam  theory  are  given 
in  [18],  [20],  and  [21].  For  the  sake  of  completeness,  we  briefly  describe  the  topics  of  the  theory  that 
we  utilize  in  this  work. 

The  Hermite-Gaussian  modes  comprise  a  family  of  solutions  to  the  two-dimensional  paraxial 
wave  equation  in  free  space.  The  Hermite-Gaussian  mode  4>n(x,z)  of  order  n  is  of  the  form 
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In  Equation  (102),  Hn  is  the  Hermite  polynomial  of  order  n.  The  quantity  w(z)  is  referred  to 
as  the  spot  size  of  the  beam,  and  is  a  measure  of  the  width  of  a  Gaussian  beam.  The  expression  for 
w(z)  is  of  the  form 


’  (103) 

where  wq  is  the  minimum  spot  size,  and  A  is  the  wavelength  of  the  light.  The  quantity  q(z)  in  (102) 
is  referred  to  as  the  complex  beam  parameter,  and  its  evolution  is  given  by: 

,  .  inwk  fir\A\ 

q(z)  =  z  -f  =  z  +  q0  .  (!04) 

As  can  be  seen  from  (102),  (103),  and  (104),  the  Hermite-Gaussian  modes  are  completely  determined 
once  values  for  wo  and  A  are  specified. 

The  Hermite-Gaussian  modes  given  by  (102)  have  several  useful  properties.  The  modes  0n(x,z) 
obey  an  orthonormality  relation 
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where  Snm  is  the  Kronecker  delta.  Moreover,  the  Hermite-Gaussian  modes  provide  a  complete  set 
of  basis  functions  characterized  by  the  single  complex  parameter  wq  [19].  As  a  result,  they  can  be 
used  to  expand  a  solution  u{x,  y,  z)  of  the  free  space  paraxial  wave  equation  as 


u(x,y,z)  =  y,  Y.  cnm(z)(j>n(x,  z)<i>m(y,  z). 

n  m 


(105) 
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Multiplying  both  sides  of  (105)  by  4>*n(x,z)4>*m(y,z)  and  integrating  over  the  transverse  plane 
yields 


/oo  poo 

/  u(x,y,z)<j>*n(x,z)^>*m(y,z)dxdy.  (106) 

-oo  J—  OO 

Thus,  a  solution  u{x^y^z)  of  the  paraxial  wave  equation,  characterized  by  u>o,  can  be  expanded  as 
in  (105)  where  the  coefficients  used  in  the  expansion  are  given  by  (106).  If  u(x,y,z)  is  a  solution  to 
the  free-space  paraxial  wave  equation,  then  the  coefficients  cnm(z)  will  be  indepenedent  of  z. 

The  spherical  mirror  transforms  the  complex  beam  parameter  q(z)  used  in  the  construction  of 
the  Hermite-Gaussian  modes.  The  transformation  law  for  a  spherical  mirror  is  given  by 

-  =  (107) 

42  4i  r 

where  r  is  the  mirror  radius.  In  (107),  the  quantity  41  is  the  value  of  q(z)  of  the  wave  incident  on 
the  curved  mirror,  while  the  value  of  the  complex  beam  parameter  for  the  wave  reflecting  from  the 
curved  mirror  is  given  by  42 . 

The  influence  of  the  spherical  mirror  can  be  obtained  by  expanding  a  solution  of  the  paraxial  wave 
equation  in  terms  of  the  Hermite-Gaussian  modes  at  the  curved  mirror,  transforming  the  modes, 
and  then  reconstructing  the  wave  consisting  of  the  transformed  modes.  In  the  following  section, 
we  utilize  the  properties  of  Gaussian  beams  to  formulate  a  continuum  model  describing  the  electric 
field  inside  the  optical  cavity  of  the  COIL.  More  generally,  the  effect  of  other  optical  elements  on 
Hermite-Gaussian  beams  can  be  described  by  so-called  ABCD  transformation  rules  [18,  see  Chapter 
20]  applied  to  the  complex  beam  parameter,  q. 

10-2  Optical  Cavity  Model 

We  are  interested  in  laser  intensity  after  one  round-trip  through  the  optical  cavity.  Thus,  we  specify 
an  initial  electromagnetic  wave  profile  at  z  =  0,  moving  upwards  from  the  outcoupling  mirror  (Mi, 
see  Figure  33).  The  solution  is  propagated  through  increasing  2  until  the  wave  reaches  the  spherical 
mirror,  M2.  At  the  spherical  mirror,  the  propagating  wave  is  reflected,  and  the  wave  propagates 
towards  decreasing  values  of  z  until  it  arrives  at  the  horizontal  mirror  again.  Upon  arriving  at  the 
outcoupling  mirror,  a  portion  of  the  wave  energy  is  extracted  due  to  energy  transmission  through 
the  mirror.  In  the  RADICL  device  (a  research  and  development  version  of  COIL,  [17,  see  Chapter 
1]  )  approximately  ten  percent  of  the  wave  energy  is  transmitted  through  the  outcoupling  mirror, 
while  roughly  ninety  percent  of  the  energy  is  reflected  back  into  the  optical  cavity. 

We  are  now  in  a  position  to  describe  a  continuum  model  for  the  electric  field  inside  the  optical 
cavity.  Since  we  start  with  an  initial  condition  at  2?  =  0  and  propagate  the  solution  toward  the 
spherical  mirror,  the  first  part  of  our  model  describes  wave  propagation  toward  increasing  values  of 
z . 


10.2.1  Propagation  Toward  Increasing  z 

The  first  part  of  our  model  describes  wave  propagation  from  M\  to  M2,  toward  increasing  values  of 
z.  As  a  result,  wave  behavior  is  modeled  by 


du 

~Fz 


f  d2u  92u\ 
\<9z2  +  dy2  ) 


+  [a  —  ik(n  -  1  )]u, 


(k  >  0)  . 


(108) 


An  initial  condition  is  specified  at  the  horizontal  mirror.  Care  must  be  taken  in  specifying  the 
initial  condition.  Our  aim  is  to  utilize  the  properties  of  Gaussian-beam  reflection  at  the  spherical 
mirror.  Recall  that  the  Hermite-Gaussian  basis  functions  presented  in  the  previous  section  depend 
on  the  beam  spot  size  w(z).  In  this  work,  we  specify  that  the  minimum  spot  size  u>0  is  located  at 
the  outcoupling  mirror  initially.  We  use  what  is  known  as  the  99%  rule  to  initially  specify  Wq  at 
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the  flat  mirror  M\ .  This  rule  is  used  as  a  design  criterion  in  determining  proper  aperture  size.  An 
aperture  of  half-length  ar  will  pass  just  over  99%  of  the  energy  of  a  Gaussian  beam  if  ar  satisfies 


CLf  — 


(109) 


where  w  is  the  spot  size  at  the  aperture. 

The  relationship  given  in  (109)  is  the  99%  rule.  We  utilize  (109)  by  specifying  the  half-length  of 
the  outcoupling  mirror  Mi,  and  then  determine  w0  according  to 


wo 


2  Mr 

7T 


(110) 


where  Mr  is  the  half-length  of  the  flat  mirror.  We  then  specify  an  initial  condition  of  the  form 


u{x,y,  0)  =u0(x,y,w0), 


(111) 


where  wq  is  the  spot  size  resulting  from  the  99%  rule.  In  the  event  that  the  half-lengths  of  the  flat 
mirror  are  unequal  in  the  x  and  y  directions,  we  utilize  (110)  and  specify  the  spot  sizes  in  the  x  and 
y  directions  separately. 

While,  in  principle  the  wave  diffuses  to  great  distances  in  the  transverse  directions,  in  practice 
the  energy  of  the  propagating  wave  is  confined  near  the  optical  axis.  Thus,  we  specify  a  to  be 
several  multiples  of  the  half-length  of  Mi  in  the  x  direction,  and  b  to  be  several  multiples  of  the 
half-length  of  Mi  in  the  y  direction.  The  outcoupling  mirror  is  then  located  in  the  interior  of  the 
region  [-a,  a]  x  [ -b,b ].  Choosing  a  and  b  sufficiently  large  allows  us  to  specify  that  the  electric 
field  remain  at  zero  along  the  transverse  boundary.  In  addition,  due  to  the  geometry  of  the  COIL 
apparatus,  we  assume  that  the  electric  field  is  symmetric  about  y  =  0.  The  resulting  boundary 
conditions  we  prescribe  are  of  the  form 


u(—a,  y,  z )  =  u(a,  y,  z )  =  u{x ,  6,  z)  =  0  , 


and 


du(a;,0,z)  _ 
dy 


(112) 

(113) 


The  partial  differential  equation  model  above  is  applied  for  wave  propagation  from  the  horizontal 
mirror  to  z  —  D. 


10.2.2  Mirror  Transformation 

Upon  arriving  at  z  =  D,  the  propagating  wave  encounters  the  spherical  mirror.  Recall  that  a  spher¬ 
ical  mirror  transforms  the  complex  beam  parameter  q(z)  used  in  the  construction  of  the  Hermite- 
Gaussian  modes.  In  order  to  incorporate  the  appropriate  transformation  into  our  model,  we  ex¬ 
pand  the  solution  u(x,y,D)  in  terms  of  the  Hermite-Gaussian  basis  functions.  That  is,  we  expand 
u{x,y,D)  as 

^(®»Vi^)  =  ^  ^  ^  )  Cnm^n(-E)  i^)<ftm(l/;  -P)i  (H4) 

n  m 

where  <j>n  and  <j>m  are  the  Hermite-Gaussian  modes  given  by  (102,  103,  104).  Recall  that  the  coeffi¬ 
cients  used  in  the  expansion  are  computed  according  to  (106),  viz.  : 


Cnm  — 


u(x,  y,  D)(j>*n{x,  D)4>*m(y,  D)dxdy. 


If  we  denote  the  value  of  the  complex  beam  parameter  at  z  =  D  by  q{D)  =  qi  ,  then,  by  utilizing 
(107),  the  transformed  beam  parameter  92  due  to  reflection  is  given  by 


A  - 1  _  ? 

92  qi  T  ' 


(115) 
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where  r  is  the  radius  of  the  spherical  mirror,  M2.  The  transformed  Hermite- Gaussian  modes  are  of 
the  form 


0>n(x,z) 


. kx 2 


(116) 


Recombining  the  transformed  Hermite-Gaussian  modes  yields  the  transformed  solution  u(x,  y,  D) 
due  to  reflection.  The  resulting  expression  for  u(x,y,  D)  is  of  the  form 


u(x,y,D)  =  EE  Cnm  4>n(x,D)4>m(y,D). 

n  m 


(117) 


10.2.3  Propagation  Toward  Decreasing  z 

After  reflection  at  z  =  D,  the  transformed  wave  then  propagates  toward  decreasing  values  of  z . 
Thus,  wave  propagation  after  reflection  is  described  by  the  downward  paraxial  wave  equation  of  the 
form 


I'2ffi(^  +  l?)+[“<X’!'’2)  +  ’l‘’l("_1)l9’  {k<°h  (U8) 

The  transformed  solution  at  z  =  D  is  used  as  the  initial  data  for  (118).  That  is,  we  specify  that 

g(x,y,D)  =  u(x,y,D).  (119) 

As  in  the  case  of  upward  propagation,  we  require  that  the  electric  field  remain  at  zero  along  the 
transverse  boundary.  Moreover,  we  assume  that  the  electric  field  is  symmetric  about  y  —  0  and 
impose  the  boundary  conditions  (112),  (113). 


10.2.4  Energy  Extraction 

After  reflecting  from  the  spherical  mirror,  the  wave  propagates  through  decreasing  z  values  until 
it  arrives  back  at  the  outcoupling  mirror  M\.  The  modulus  of  the  electric  field  incident  on  the 
outcoupling  mirror  is  given  by  \g{x,y,  0)|.  In  order  to  completely  describe  wave  behavior  for  one 
round-trip  through  the  optical  cavity,  we  now  model  the  extraction/ reflection  of  wave  energy  that 
occurs  at  the  outcoupling  mirror.  As  noted  above,  in  the  RADICL  device  approximately  ten  percent 
of  the  wave  energy  is  transmitted  through  the  outcoupling  mirror,  so  that  roughly  ninety  percent  of 
the  wave  energy  is  reflected  from  the  outcoupling  mirror  back  into  the  optical  cavity. 

Recall  that  wave  reflection  from  a  spherical  mirror  transforms  the  complex  beam  parameter  q{z) 
used  in  the  specification  of  the  Hermite-Gaussian  modes.  As  can  be  seen  from  (107),  the  parameter 
q(z)  is  unchanged  in  the  case  of  reflection  from  a  mirror  with  infinite  radius  of  curvature.  As  a  result, 
to  determine  the  wave  profile  after  one  round-trip  through  the  optical  cavity,  we  simply  decompose 
g(x)  y,  0)  into  two  parts.  The  first  part  is  that  portion  of  g(x,  y,  0)  that  is  incident  on  the  outcoupling 
mirror.  The  second  part  is  the  portion  of  g(x,y,0)  that  misses  the  flat  mirror. 

To  incorporate  the  transmission  of  energy  through  the  outcoupling  mirror,  define  the  energy 
Su(z)  associated  with  a  wave  profile  u(x,y,z)  according  to 


rb  ra 

£u(z)  =  /  /  \u(x,y,z)\2dxdy. 

J~b J-a 


Define  the  function  gM(x,y> 0)  by 


9m{x,v,  0) 


g(x,y,0)  if  (x,y)  €  P(Mi); 
0,  otherwise, 


(120) 
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where  V(Mi)  is  the  subset  of  the  region  [-a, a]  x  [0,6]  occupied  by  the  flat  mirror.  In  a  similar 
fashion,  define  gii(xyy,  0)  according  to 


9r(x,v,0) 


g{x,y,0)  if  (x,y)  ^  P(Mj); 
0,  otherwise. 


For  the  sake  of  simplicity,  we  assume  that  energy  not  incident  on  the  outcoupling  mirror  is 
absorbed  by  the  optical  cavity  wall.  As  a  result,  the  energy  associated  with  gR.(x,y)  is  completely 
absorbed,  and  none  of  its  energy  is  reflected  back  into  the  optical  cavity. 

The  requirement  that  ten  percent  of  the  energy  of  the  wave  profile  incident  on  M\  is  transmitted 
through  the  mirror  results  in  the  specification  that 

/oo  rOO  r  CO  r  OO 

/  \uR{x,y,0)\2dxdy  =  0.9  /  /  \gM(x,y,0)\2dxdy, 

-oo  J — oo  J  — oo  J  —oo 


where  ur(x,  y,  0)  is  the  wave  profile  after  one  round-trip  through  the  optical  cavity.  The  90%  reflec¬ 
tion  requirement  is  applied  pointwise,  so  that  the  round-trip  solution  u#(:r,y,  0)  is  found  according 
to 


un(x,y,0)  =  x/09  9m(x,v, 0). 


10.3  Numerical  Implementation 

In  order  to  obtain  numerical  solutions  to  the  initial  boundary  value  problem  presented  in  the  previous 
section,  we  construct  a  uniform  finite  difference  grid  on  the  physical  domain  represented  in  Figure 
33.  The  step-size  in  the  z  direction  is  denoted  by  A z.  The  step-size  in  the  x  and  y  directions  is 
denoted  by  h.  Evaluation  points  in  the  x  direction  are  in  the  interval  [-a,  a].  Similarly,  y  values  are 
in  the  interval  [0,6]. 


10.3.1  Propagation  Toward  Increasing  z 


As  we  are  first  interested  in  propagation  toward  increasing  values  of  z,  our  finite  difference  scheme 
for  upward  propagation  must  approximate  the  initial  boundary  value  problem  given  by  (108),  (111)- 
(113).  To  this  end,  the  first  order  derivative  in  the  z  direction  is  approximated  as 


d  ,  N  u(v,z  +  Az)-u(v^) 
&“<•’  z)  “ - Tz - 


(121) 


The  second  order  derivatives  in  the  transverse  directions  are  approximated  by 


g^u(x,.,.)* - -2 - ,  and 

d2  ,  s  u(-,y-h,-)  -2u(-,y,  •)  +  «(•,  y  +  V) 
y,  •)  » - ^ - • 


(122) 

(123) 


In  order  to  simplify  the  resulting  finite  difference  approximation,  we  introduce  the  quantities  v 
and  cfj  given  by 


tAz 


v 


and 


2kh2 ' 

clj  =Az«j  -  *  k  ( ni,j  -  !))• 


(124) 

(125) 


The  resulting  finite  difference  scheme  we  utilize  at  each  evaluation  point  (x{,yj,zn)  is  a  Crank- 
Nicholson  scheme  of  the  form 


-  K+ 1  j  +  (2  +  4i/  +  c",) 


(126) 
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Convergence  properties  of  (126)  are  given  in  [14]. 

The  symmetry  condition  given  by  (113)  is  incorporated  into  (126)  in  the  case  of  j  —  0.  That  is, 
for  j  =  0,  the  scheme  given  by  (126)  becomes 


[vu 


71+1 

i — 1,0 


+K?i1,o  +  (2  -  4*  -  c#1)^1]  +  2[<]  = 

[— —  vu7+i,  o  +  (2  +  4*  +  <oKo]  —  2[^'wil]. 


(127) 


Iterations  for  upward  propagation  are  done  using  (126)  and  (127).  The  iterations  are  started 
using  the  initial  condition  uo(xi, ^,0). 


10.3.2  Mirror  Transformation  Implementation 

Utilizing  the  Crank-Nicholson  scheme  of  the  previous  section,  one  obtains  a  numerical  solution  to  the 
initial  boundary  value  problem  describing  upward  wave  propagation.  At  an  axial  distance  of  z  —  D^ 
the  propagating  wave  encounters  the  spherical  mirror  M2.  Thus,  we  now  transform  the  solution 
at  z  =  D  in  order  to  incorporate  the  properties  of  Gaussian  beam  reflection  into  the  numerical 
approximation.  That  is,  we  expand  the  solution  u{xi,yj,D )  in  terms  of  the  Hermite-Gaussian 
modes  as 

Nm  Nm 

W(^i>  Vj )  -O)  =  EE  Cnm$n{xii  &)4>rn{yj )  D))  (128) 

n=0  771—0 

where  Nm  is  the  order  of  the  highest  order  mode  used  in  the  expansion. 

In  order  to  utilize  the  orthonormality  properties  of  the  Hermite-Gaussian  modes,  we  reflect 
the  numerical  solution  about  y  =  0  to  get  a  numerical  solution  with  support  on  the  rectangle 
[-a,  a]  x  [-6,6].  The  expansion  coefficients  cnm  are  then  found  according  to 

Cnm  =  / 6  P  u(x,y,D)<t>*n(x,D)4>*m(y,D)dxdy.  (129) 

J-b J-a 

The  integral  in  (129)  is  approximated  via  a  double  Riemann  sum. 

Having  obtained  the  coefficients  cnm  needed  in  the  expansion,  the  Hermite-Gaussian  modes  are 
transformed  according  to  Equation  (116).  Recombining  the  transformed  Hermite-Gaussian  modes 
yields  a  transformed  numerical  solution  u(xi,yj,D)  of  the  form 

Nm  Nm 

u(xi,yj,D)  =  EE  C  71 771 0  7i  {xuD)4m{yj,D\  (130) 

71—0  771=0 

where  <^n  is  the  transformed  mode  of  order  n. 


10.3.3  Propagation  Toward  Decreasing  z 


We  now  shift  our  attention  to  the  propagation  of  the  reflected  wave  back  to  the  outcoupling  mirror. 
The  numerical  scheme  needed  for  downward  propagation  is  easily  obtained  by  reversing  the  sign  of 
k ,  and  results  in  a  Crank-Nicholson  scheme  of  the  form 


[-■*££/+  -  &  +  (2+4*- <+>  K+1]  -  KA  ]  - 


+  ”9?+  i,j  +  (2  -  4 u  +  c^glj]  +  +  [vg?ij+1]. 


(131) 


Convergence  properties  of  (131)  are  discussed  in  [14]. 

As  before,  the  symmetry  condition  specified  by  (113)  is  incorporated  into  (131)  for  the  case  of 
j  =  0.  That  is,  in  the  case  of  j  =  0,  (131)  becomes 


(132) 


The  remaining  iterations  are  done  with  (131)  and  (132).  The  transformed  numerical  solution 
u(xi,2/j,D)  due  to  reflection  is  used  as  initial  data. 
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10.3.4  Energy  Extraction  Implementation 

In  order  to  obtain  the  wave  profile  after  one  round-trip  through  the  optical  cavity,  we  decompose 
g(xiyyj, 0)  into  two  parts,  as  described  in  Section  10.2.4.  The  round-trip  solution  uR(xi,yj, 0)  is 
then  found  according  to 

UR{xi,yj, 0)  =  Vo$  gM(xuyj,0),  (133) 

where  0Af(®i>yj>  0)  is  that  portion  of  g{xiyy^  0)  incident  on  the  flat  mirror. 


10.3.5  Average  Two-Way  Intensity 


Source  terms  in  the  reacting  flow  model  require  an  intensity  field  in  the  optical  cavity.  Since  the 
flow  model  is  based  on  a  unit-cell  approximation,  it  is  appropriate  to  supply  an  intensity  field  that 
is  averaged  (in  the  z  direction) .  For  this  purpose  we  select  a  number  of  ^-planes  and  compute  the 
average  two-way  intensity  field  as: 


/avg  (%i  >  Vj )  —  jy 


(134) 


where  {k(l)  |  l  =  1,.. .  ,NP}  is  the  index  set  labelling  the  z  planes,  and  u  ( g )  is  the  upward  ( resp . 
downward)  electric  field  solution. 


10.4  The  Coupled  Model 

As  previously  noted,  analyses  codes  for  continuous  wave  chemical  lasers  can  be  viewed  as  coupled 
fluid  mechanics  and  laser  power  extraction  codes.  In  the  present  study  we  have  used  AeroSoft’s 
QASVvA  code  to  solve  for  the  flow  field  (u/),  given  the  average  two-way  intensity  ( ui )  in  the  optical 
cavity.  The  calculations  described  in  Section  10.3  were  implemented  in  a  C++  code  so  that  an 
initial  electric  field  at  the  flat  outcoupling  mirror  was  propagated  one  round-trip  through  the  optical 
cavity.  As  part  of  this  calculation  we  store  a  new  average  two-way  intensity  in  the  optical  cavity 
[see  Equation  (134)  ].  With  a  slight  abuse  of  notation  we  write  the  system  abstractly  as: 

ui-'Ui{uf,ui)  =  0,  (135) 

Kfiuf.Ui)  =  0  (136) 

In  the  COIL  application  the  variable  it*  is  interpreted  as  as  pair  of  fields:  the  initial  electric  field 
at  the  outcoupling  mirror,  and  the  average  two-way  intensity  in  the  active  gain  region.  The  residual 
for  the  chemically  reacting  flow  (136)  depends  on  the  latter  of  these  fields. 


10.4.1  Model  Problem 

In  order  to  study  solution  procedures  for  coupled  systems  of  the  form  (135,  136)  we  constructed  a 
model  problem  based  on  one-dimensional  gas  dynamics  with  heat  addition  [12].  This  model  problem 
was  also  used  to  study  sensitivity  calculations  for  such  coupled  systems. 

10.4.2  Numerical  Results 

Two  numerical  implementations  of  the  coupled  COIL  system  (135,  136)  were  studied.  At  AeroSoft 
Inc.  calculations  were  based  on  the  QASVvA  reacting  flow  code  coupled  to  a  ray-tracing  optics 
model  described  by  Crowell  [2]. 

The  numerical  studies  at  ICAM  were  based  on  GASVvA  loosely  coupled  to  our  C++  implementa¬ 
tion  of  the  paraxial  optics  model.  At  the  end  of  the  effort  there  were  unresolved  issues  in  the  coupled 
codes  so  that  no  converged  results  are  available.  Figure  34  shows  a  typical  twoway  intensity  field 
after  several  iterations.  It  appears  that  the  beam  parameters  used  to  define  the  Her  mite- Gaussian 
beam  family  are  not  appropriate  and  attempts  to  tune  these  parameters  have  not  yet  been  successful. 
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Figure  34:  Two  Way  Intensity  Field 

10.5  Semigroup  Formulation 

From  a  mathematical  point  of  view  the  partial  differential  equation  system  given  by  Equations  (108, 
111,  112,  113)  is  a  formal  model  since  it  is  not  precisely  stated  what  regularity  an  initial  function 
must  have  or  what  regularity  the  solution  will  have.  For  example,  given  the  solution  at  station  z 
does  the  energy  expression  (120)  have  a  finite  value?  The  analysis  in  Camphouse’s  PhD  thesis  [14] 
provides  a  framework  and  answers  for  these  questions.  Here  we  reproduce  some  of  the  material  from 
Camphouse’s  work  [14,  see  Chapter  3,  26  -  31.].  The  presentation  is  slightly  altered  in  a  way  that 
suggests  an  alternative  numerical  treatment. 

Consider  the  rectangular  domain  $4  =  [-a,  a]  x  [—6  x  b ]  and  let  %  =  £2(^)1  the  square  integrable 
functions  on  the  specified  domain.  Consider  the  upward  wave  model  (108)  with  a  =  0  and  n  =  1 
and  define  the  operator  Aq  :  V(Aq)  C  H  H  by 

1  ( d2u  d2u\  f  . 

AoU  =  ~2k{w  +  WJ  •  (  37) 

From  [14,  Theorem  3.1.1]  we  have  The  operator  ^4o  defined  by  (137)  generates  a  Co-group  of 
unitary  operators  {<So(^)|  —  00  <  z  <  00}  on  the  space  H.  In  simple  terms,  for  some  fixed  z ,  the 

group  element  *So(£)  is  the  solution  operator  that  takes  an  initial  electric  field  u0(x,  2/)  at  2  =  0  and 

maps  it  to  the  (unique)  solution  to  the  boundary  value  problem  (108,  111-113)  at  z  =  z.  The  utility 
of  the  semigroup  formulation  is  that  it  establishes  that  there  is  a  unique  solution  to  the  boundary 
value  problem  and  that  the  solution  depends  continuously  in  the  initial  data.  Note  that  Aq  models 
the  free-space  behavior  of  the  paraxial  waves. 

If  the  gain  is  non-zero  or  the  index  of  refraction  is  not  unity  then  there  is  an  additional  term  in 
the  paraxial  wave  equation  (108).  For  this  purpose  we  require  that  the  complex- valued  function 

g{x,y)  =  a(x,y)  -  ik{n(x,y)  -  1) 

be  in  W  and  define  the  operator  A\  on  V(A\)  =  T>(Ao)  C  W,  W  by 

Axu  =  Aou  +  5U,  (138) 

where,  as  expected,  g  u  means  pointwise  multiplication.  From  [14,  Theorem  3.1.2]  we  have  The 
operator  A{  defined  by  (138)  generates  a  Co-semigroup  of  operators  {<Si(*)|  0  <  z  <  00}  on 
Moreover,  the  semigroup  action  is  given  by 

Si(z)u  =  exp(<7  2)iSo(z)u  •  (139) 

The  representation  (139)  means  that  we  can  find  a  solution  to  the  problem  with  gain  by  first  find¬ 
ing  the  solution  with  no  gain  and  multiplying  the  resulting  field  pointwise  by  the  value  exp  [g(x,  y)  z). 
This  idea  can  be  used  as  a  basis  for  an  alternative  numerical  formulation.  The  calculations  we  have 
in  mind  here  are  very  similar  to  those  reported  in  [19]. 
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11  Personnel  Supported 

Key  personnel  during  the  Phase-I  effort  are  Dr.  William  Eppard,  Dr.  Andrew  Godfrey  and  Dr.  Michael 
Applebaum  at  Aerosoft  and  Dr.  Eugene  Cliff,  Dr.  John  Burns  and  Dr.  Chris  Camphouse  with  ICAM 
at  Virginia  Tech. 


12  Publications 

Reference  [4]  and  [3]  above. 

13  Interactions/Transitions 

13.1  Participation/Presentations  at  Meetings 

Presentations  at  the  31st  Plasmadynamics  and  Lasers  Conference  in  Denver  (19-22  of  June  2000) 
and  the  33rd  AIAA  Plasmadynamics  and  Lasers  Conference  in  Maui  (20-23  of  May). 

13.2  Transitions 

The  improved  modeling  capabilities  in  GASP  will  provide  an  immediate  impact  on  COIL  analysis 
efforts  at  AFRL/DE.  TRW  and  the  Schafer  Corp.  are  also  using  GASP  for  chemical  laser  analysis. 


14  New  Discoveries,  Inventions,  or  Patent  Disclosures 

None. 

1 5  Honors  /  Awards 

None. 


61 


References 


[1]  R.  C.  Buggeln,  S.  Shamroth,  A.  I.  Lampson,  and  P.  G.  Crowell.  “Three-Dimensional  (3-D) 
Navier-Stokes  Analysis  of  the  Mixing  and  Power  Extraction  in  a  Supersonic  Chemical  Oxygen 
Iodine  Laser  (COIL)  With  Transverse  I2  Injection”.  AIAA  Paper  94-2435, 25th  Plasmadynamics 
and  Lasers  Conference,  July,  1993. 

[2]  P.  G.  Crowell.  “A  Stable  Resonator  Geometric  Optics  Model  for  Gain  Saturation  and  Power 
Extraction  in  COIL  Devices”.  Contract  F29601-93-C-0028,  CDRL  A005,  Phillips  Laboratory, 
Februrary,  1995. 

[3]  W.  M.  Eppard,  W.  D.  McGrory,  and  M.  P.  Applebaum.  “The  Effect  of  Water- Vapor  Con¬ 
densation  and  Surface  Catalysis  on  COIL  Performance”.  AIAA  Paper  2002-2132,  33rd  AIAA 
Plasmadyanmics  and  Lasers  Conference,  May  20-23,  2002. 

[4]  W.  M.  Eppard,  W.  D.  McGrory,  A.  G.  Godfrey,  E.  M.  Cliff,  and  J.  T.  Borggaard.  “Recent 
Advances  in  Numerical  Techniques  for  the  Design  and  Analysis  of  COIL  Systems” .  AIAA  Paper 
2000-2576,  31s<  AIAA  Plasmadyanmics  and  Lasers  Conference,  June  19-22,  2000. 

[5]  A.  G.  Godfrey  and  E.  M.  Cliff.  “Direct  Calculation  of  Aerodynamic  Force  Derivatives:  A 
Sensitivity  Equation  Approach”.  AIAA  Paper  98-0393,  36t/l  AIAA  Aerospace  Sciences  Meeting 
and  Exhibit,  January  12-15, 1998. 

[6]  A.  G.  Godfrey,  W.  M.  Eppard,  and  E.  M.  Cliff.  “Using  Sensitivity  Equations  for  Chemically 
Reacting  Flows”.  AIAA  Paper  98-4805,  7th  AIAA/USAF/NASA/ISSMO  Symposium  on  Mul¬ 
tidisciplinary  Analysis  and  Optimization,  September  2-4,  1998. 

[7]  J.  O.  Hirschfelder,  C.  F.  Curtiss,  and  R.  B.  Bird.  Molecular  Theory  of  Gases  and  Liquid. 
John  Wiley  and  Sons,  Library  of  Congress  CCN  54-7621,  1954. 

[8]  M.  Hishida,  N.  Azami,  K.  Iwamoto,  W.  Masuda,  H.  Fujii,  T.  Atsuta,  and  M.  Mura.  “Flow  and 
Optical  Fields  in  a  Supersonic  Flow  Chemical  Oxygen-Iodine  Laser”.  AIAA  Paper  97-2391, 
28th  Plasmadynamics  and  Lasers  Conference,  June  23-25  1997. 

[9]  T.  J.  Madden.  “Personal  Communication”,  January  2001. 

[10]  W.  Masuda,  M.  Satoh,  and  H.  Yamada.  “Effects  of  Water  Vapor  Condensation  on  Performance 
of  Supersonic  Flow  Chemical  Oxygen-Iodine  Laser”.  JSME  Int.  J.,  39(2):273,  1996. 

[11]  J.  T.  Borggaard,  E.M.  Cliff  and  W.  M.  Eppard.  “Working  Notes  on  Sensitivity  Modifications 
to  stable. F  ”  ICAM  Report  01-0 f -00.  VPI  &  SU,  Blacksburg,  VA,  April,  2000. 

[12]  J.T.  Borggaard  and  E.M.  Cliff.  “Sensitivity  Analysis  for  Chemical  Laser  Design:  A  Model 
Problem.”  ICAM  Report  01-06-00.  VPI  &  SU,  Blacksburg,  VA,  June,  2000. 

[13]  R.  C.  Buggeln  and  S.  Shamroth  and  A.  I.  Lampson  and  P.  G.  Crowell.  “Three-Dimensional 
(3-D)  Navier-Stokes  Analysis  of  the  Mixing  and  Power  Extraction  in  a  Supersonic  Chemical 
Oxygen  Iodine  Laser  (COIL)  With  Transverse  I2  Injection.”  AIAA  Paper  94-2435.  June  1994. 

[14]  R.  Chris  Camphouse.  Modeling  and  Numerical  Approximations  of  Optical  Activity  in  the  Chem¬ 
ical  Oxygen-Iodine  Laser.  PhD  Thesis,  Department  of  Mathematics,  VPI  &  SU,  Blacksburg, 
VA,  July  2001. 

[15]  P.G.  Crowell.  A  Stable  Resonator  Geometric  Optics  Model  for  Gain  Saturation  and  Power 
Extraction  in  COIL  Devices,  Phillips  Laboratory,  Contract  F29601-93-C-0028,  CDRL  A005. 

[16]  Andrew  G.  Godfrey,  William  M.  Eppard,  Eugene  M.  Cliff.  “Using  Sensitivity  Equations  for 
Chemically  Reacting  Flows.”  AIAA  Paper  98-4805  .  7th  AIAA/USAF/NASA/ISSMO  Sym¬ 
posium  on  Multidisciplinary  Analysis  and  Optimization.  St.  Louis,  Missouri.  2-4  September, 
1998. 


62 


[17]  Timothy  J.  Madden.  Computational  Fluid  Dynamics  Methodologies  for  Simulation  of  Chemical 
Oxygen-Iodine  Laser  Flowfields.  PhD  Thesis,  Department  of  Aeronautical  and  Astronautical 
Engineering,  University  of  Illinois  at  Urbana-Champaign,  1997. 

[18]  Anthony  E.  Siegman.  Lasers.  University  Science  Books,  Mill  Valley,  CA,  1986. 

[19]  Anthony  E.  Siegman,  Edward  A.  Sziklas.  “Mode  Calculations  in  Unstable  Resonators  with 
Flowing  Saturable  Gain.  1:  Hermite-Gaussian  Expansion.”  Applied  Optics.  Vol.  13,  No.  12. 
December,  1974. 

[20]  Joseph  T.  Verdeyen.  Laser  Electronics ,  Third  Edition,  Prentice  Hall,  1995. 

[21]  Amnon  Yariv.  Introduction  to  Optical  Electronics.  Holt,  Rinehart  and  Winston,  1976. 


63 


